# a script that runs the structural decomposition analysis (SDA) of the changes in GDP, industry value added and employment induced by bilateral trade creation, diversion and contraction
# method of the bilateral allocation of market share losses to market share gains: simplified, bi-proportional
# version: Dec 2019
# ENSURE that Octave can find the script, use the "addpath" command: addpath ("/[insert the path to your Octave folder]/Octave")

# enter year 0 (2005 <= y0 <= 2014) and year 1 (2006 <= y1 <= 2015)
year = input ("enter year 0: ")
y0 = num2str (year);
year = input ("enter year 1: ")
y1 = num2str (year);
# construct an aggregation matrix S0 to aggregate the data for Mexico (MX1 - Mexico - Activities excluding Global Manufacturing, MX2 - Mexico - Global Manufacturing activities) and China (CN1 - China - Activities excluding export processing, CN2 - China - Export processing activities) into the rows/columns for MEX and CHN, reducing the number of country categories from 69 to 65
# S0 is a 2484 x 2340 sparse matrix
S0 = sparse (2484, 2340);
S0(1:792, 1:792) = speye(792);
S0(2341:2412, 793:828) = [speye(36); speye(36)];
S0(829:1476, 829:1476) = speye(648);
S0(2413:2484, 1477:1512) = [speye(36); speye(36)];
S0(1513:2340, 1513:2340) = speye(828);
# construct an aggregation matrix Sf to aggregate final demand categories within each country
# Sf is a 390 x 65 sparse matrix
Sf = kron(speye(65),ones(6,1));
# load and format ICIO for year 1
filename1 = ["/[insert the path to the OECD ICIO Octave files]/ICIO2018_" y1 "_octave"];
load (filename1, "Z0", "F0", "TLSi0", "TLSf0", "v0", "x0")
Z1 = S0'*Z0*S0;
F1 = S0'*F0*Sf;
TLSi1 = TLSi0*S0;
TLSf1 = TLSf0*Sf;
x1 = x0*S0;
v1 = v0*S0;
clear Z0 F0 TLSi0 TLSf0 v0 x0
# load and format ICIO for year 0
filename0 = ["/[insert the path to the OECD ICIO Octave files]/ICIO2018_" y0 "_octave"];
load (filename0, "Z0", "F0", "TLSi0", "TLSf0", "v0", "x0")
Z0 = S0'*Z0*S0;
F0 = S0'*F0*Sf;
TLSi0 = TLSi0*S0;
TLSf0 = TLSf0*Sf;
x0 = x0*S0;
v0 = v0*S0;
clear S0 Sf
# load employment data for years 0 and 1
filename = ["/[insert the path to the Octave files with the processed OECD TiM data]/EMPN_octave"];
e = load (filename, ["e_" y0],["e_" y1]);
# e is a structure array with two fields
e0 = e.(["e_" y0]);
e1 = e.(["e_" y1]);
clear e

# compute the matrix of technical coefficients A for year 0
A0 = Z0./x0;
# matrix A may contain "NaN" values that usually result from division of zero by zero; this is typical for the columns that correspond to "Private households with employed persons" (D97T98) if this industry is not reported to produce any output
# replace "NaN" in A with zeros
A0(isnan(A0)==1) = 0;
# compute the Leontief inverse L for year 0
L0 = inv(eye(2340)-A0);
# compute A and L for year 1
A1 = Z1./x1;
A1(isnan(A1)==1) = 0;
L1 = inv(eye(2340)-A1);

# construct aggregation matrices
# the industry-wise aggregation matrix Sk (K = number of countries)
Sk = kron(ones(1,65),speye(36));
# the country-wise aggregation matrix Sn (N = number of industries)
Sn = kron(speye(65),ones(36,1));

# prepare component matrices for SDA, for year 0
# compute a row vector of total domestic final demand by country (dimension: 1 x K)
f0 = sum(F0,1);
# compute a matrix of industry (product) structure of domestic final demand, by country (dimension: N x K)
Find0 = (Sk*F0)./f0;
# compute a matrix of the origin of domestic final demand, by country (dimension: KN x K)
Fcou0 = F0./(Sk'*Sk*F0);
# replace "NaN" in Fcou0 with zeros
Fcou0(isnan(Fcou0)==1) = 0;
# compute row vector a, share of total intermediate inputs in output, by country-industry (dimension: 1 x KN)
a0 = sum(A0, 1);
# compute matrix Aind, industry structure of intermediate inputs, by country-industry (dimension: N x KN)
Aind0 = (Sk*A0)./a0;
# replace "NaN" in Aind0 with zeros
Aind0(isnan(Aind0)==1) = 0;
# compute matrix Acou, country origin of intermediate inputs, by country-industry (dimension: KN x KN)
Acou0 = A0./(Sk'*Sk*A0);
# replace "NaN" in Acou0 with zeros
Acou0(isnan(Acou0)==1) = 0;
# compute matrix Ti0, implied net tax rates on intermediate inputs, country by country-industry (dimension: K x KN)
Ti0 = TLSi0./(Sn'*Z0);
# replace "NaN" in Ti0 with zeros
Ti0(isnan(Ti0)==1) = 0;
# compute matrix Tf0, implied net tax rates on final products, country by country (dimension: K x K)
Tf0 = TLSf0./(Sn'*F0);
# replace "NaN" in Tf0 with zeros
Tf0(isnan(Tf0)==1) = 0;
# compute row vector ec0, employment coefficients, by country-industry (dimension: 1 x KN)
ec0 = e0'./x0;
              
# prepare component matrices for SDA, for year 1
# compute a row vector of total domestic final demand by country (dimension: 1 x K)
f1 = sum(F1,1);
# compute a matrix of industry (product) structure of domestic final demand, by country (dimension: N x K)
Find1 = (Sk*F1)./f1;
# compute a matrix of the origin of domestic final demand, by country (dimension: KN x K)
Fcou1 = F1./(Sk'*Sk*F1);
# replace "NaN" in Fcou1 with zeros
Fcou1(isnan(Fcou1)==1) = 0;
# compute row vector a, share of total intermediate inputs in output, by country-industry (dimension: 1 x KN)
a1 = sum(A1, 1);
# compute matrix Aind, industry structure of intermediate inputs, by country-industry (dimension: N x KN)
Aind1 = (Sk*A1)./a1;
# replace "NaN" in Aind1 with zeros
Aind1(isnan(Aind1)==1) = 0;
# compute matrix Acou, country origin of intermediate inputs, by country-industry (dimension: KN x KN)
Acou1 = A1./(Sk'*Sk*A1);
# replace "NaN" in Acou1 with zeros
Acou1(isnan(Acou1)==1) = 0;
# compute matrix Ti1, implied net tax rates on intermediate inputs, country by country-industry (dimension: K x KN)
Ti1 = TLSi1./(Sn'*Z1);
# replace "NaN" in Ti1 with zeros
Ti1(isnan(Ti1)==1) = 0;
# compute matrix Tf1, implied net tax rates on final products, country by country (dimension: K x K)
Tf1 = TLSf1./(Sn'*F1);
# replace "NaN" in Tf1 with zeros
Tf1(isnan(Tf1)==1) = 0;
# compute row vector ec1, employment coefficients, by country-industry (dimension: 1 x KN)
ec1 = e1'./x1;
              
# compute a matrix of change in the origin of domestic final demand, by country (dimension: KN x K)
dFcou = Fcou1 - Fcou0;
# compute a matrix of change in the origin of intermediate inputs, by country-industry (dimension: KN x KN)
dAcou0 = Acou1 - Acou0;
# extract a matrix of technical changes from dAcou that are not related to substitution between countries; these changes occur because product of industry i appears in or disappears from the intermediate consumption of industry j in country s
# construct a matrix to identify the position of the elements in dAcou that are not related to substitution between countries
dAcou_tec0 = Sk'*round(Sk*dAcou0);
# dAcou_tec0 is a matrix of ones, minus ones and zeros where ones and minus ones indicate the position of elements not related to substitution between countries (ones signify that a product appeared in intermediate consumption of country s, and minus ones signify that it disappeared)
# adjust dAcou
dAcou_tec = dAcou0.*abs(dAcou_tec0);
dAcou = dAcou0 - dAcou_tec;
clear dAcou0 dAcou_tec0
# NOTE: dAcou_tec is very likely to contain only few non-zero elements, and therefore can be made a sparse matrix, but I do not declare it sparse because it is not necessary for calculations below
# NOTE: I do not extract the same matrix of technical changes from dFcou because it is very unlikely that a product of industry i appears in or disappears from the final demand of country s
# however, I include a warning in case there are similar changes in final demand
dFcou_tec = round(Sk*dFcou);
if (max(max(dFcou_tec)) > 0)
    printf ("Beware of other changes in Fcou that are not related to substitution between countries!\n")
endif
clear dFcou_tec

# split dAcou into matrices of negative and non-negative changes
dAcou_neg = dAcou;
dAcou_neg(dAcou_neg>0)=0;
dAcou_pos = dAcou;
dAcou_pos(dAcou_pos<0)=0;
# split dFcou into matrices of negative and non-negative changes
dFcou_neg = dFcou;
dFcou_neg(dFcou_neg>0)=0;
dFcou_pos = dFcou;
dFcou_pos(dFcou_pos<0)=0;
# clear unnecessary matrices
# clear A0 A1

# construct a (sparse) permutation matrix for the RAS balancing
P0 = sparse(65,2340);
for k = 1:65
    P0(k, 1+36*(k-1)) = 1;
endfor
P1 = P0;
for n = 2:36
    P2 = circshift(P0, [0,n-1]);
    P1 = [P1;P2];
endfor
clear P0 P2
# matrix P1 shifts rows in such a way that country-industry classification becomes industry-country

# calculate individual terms in the decompositions that do not depend on country pair - INTERMEDIATE demand
# A*** = A(Acou*,Aind*,a*)
# A111 = A1;
A011 = Acou0.*(Sk'*Aind1*diag(a1));
A100 = Acou1.*(Sk'*Aind0*diag(a0));
# A000 = A0;
A110 = Acou1.*(Sk'*Aind1*diag(a0));
A010 = Acou0.*(Sk'*Aind1*diag(a0));
A101 = Acou1.*(Sk'*Aind0*diag(a1));
A001 = Acou0.*(Sk'*Aind0*diag(a1));
# L** = L(Acou*,Aind*,a*)
# L111 = L1;
L011 = inv(eye(2340)-A011);
L100 = inv(eye(2340)-A100);
# L000 = L0;
L110 = inv(eye(2340)-A110);
L010 = inv(eye(2340)-A010);
L101 = inv(eye(2340)-A101);
L001 = inv(eye(2340)-A001);

# 32 VALUE ADDED terms for changes in the country of origin of INTERMEDIATE products
# v***** = v(Tz*,Acou*,Aind*,a*,F*)
# v11111 = v1';
v10111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A011),1)))*L011*sum(F1,2);
v11001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A100),1)))*L100*sum(F1,2);
v10001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0),1)))*L0*sum(F1,2);
v11101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A110),1)))*L110*sum(F1,2);
v10101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A010),1)))*L010*sum(F1,2);
v11011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A101),1)))*L101*sum(F1,2);
v10011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A001),1)))*L001*sum(F1,2);
v01110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1),1)))*L1*sum(F0,2);
v00110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A011),1)))*L011*sum(F0,2);
v01000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A100),1)))*L100*sum(F0,2);
# v00000 = v0';
v01100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A110),1)))*L110*sum(F0,2);
v00100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A010),1)))*L010*sum(F0,2);
v01010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A101),1)))*L101*sum(F0,2);
v00010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A001),1)))*L001*sum(F0,2);
v11110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1),1)))*L1*sum(F0,2);
v10110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A011),1)))*L011*sum(F0,2);
v11000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A100),1)))*L100*sum(F0,2);
v10000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0),1)))*L0*sum(F0,2);
v11100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A110),1)))*L110*sum(F0,2);
v10100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A010),1)))*L010*sum(F0,2);
v11010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A101),1)))*L101*sum(F0,2);
v10010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A001),1)))*L001*sum(F0,2);
v01111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1),1)))*L1*sum(F1,2);
v00111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A011),1)))*L011*sum(F1,2);
v01001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A100),1)))*L100*sum(F1,2);
v00001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0),1)))*L0*sum(F1,2);
v01101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A110),1)))*L110*sum(F1,2);
v00101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A010),1)))*L010*sum(F1,2);
v01011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A101),1)))*L101*sum(F1,2);
v00011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A001),1)))*L001*sum(F1,2);

# 32 EMPLOYMENT terms for changes in the country of origin of INTERMEDIATE products
# e***** = e(ec*,Acou*,Aind*,a*,F*)
# e11111 = e1;
e10111 = diag(ec1)*L011*sum(F1,2);
e11001 = diag(ec1)*L100*sum(F1,2);
e10001 = diag(ec1)*L0*sum(F1,2);
e11101 = diag(ec1)*L110*sum(F1,2);
e10101 = diag(ec1)*L010*sum(F1,2);
e11011 = diag(ec1)*L101*sum(F1,2);
e10011 = diag(ec1)*L001*sum(F1,2);
e01110 = diag(ec0)*L1*sum(F0,2);
e00110 = diag(ec0)*L011*sum(F0,2);
e01000 = diag(ec0)*L100*sum(F0,2);
# e00000 = e0;
e01100 = diag(ec0)*L110*sum(F0,2);
e00100 = diag(ec0)*L010*sum(F0,2);
e01010 = diag(ec0)*L101*sum(F0,2);
e00010 = diag(ec0)*L001*sum(F0,2);
e11110 = diag(ec1)*L1*sum(F0,2);
e10110 = diag(ec1)*L011*sum(F0,2);
e11000 = diag(ec1)*L100*sum(F0,2);
e10000 = diag(ec1)*L0*sum(F0,2);
e11100 = diag(ec1)*L110*sum(F0,2);
e10100 = diag(ec1)*L010*sum(F0,2);
e11010 = diag(ec1)*L101*sum(F0,2);
e10010 = diag(ec1)*L001*sum(F0,2);
e01111 = diag(ec0)*L1*sum(F1,2);
e00111 = diag(ec0)*L011*sum(F1,2);
e01001 = diag(ec0)*L100*sum(F1,2);
e00001 = diag(ec0)*L0*sum(F1,2);
e01101 = diag(ec0)*L110*sum(F1,2);
e00101 = diag(ec0)*L010*sum(F1,2);
e01011 = diag(ec0)*L101*sum(F1,2);
e00011 = diag(ec0)*L001*sum(F1,2);

# for the GDP change
# calculate taxes less subsidies on final products mf**** = mf(Tf*,F*) for Tf1, Tf0, F1, F0
mf1111 = sum(TLSf1,1)';
mf0000 = sum(TLSf0,1)';
mf1000 = sum(Tf1.*(Sn'*F0),1)';
mf0111 = sum(Tf0.*(Sn'*F1),1)';
# 32 GDP terms for changes in the country of origin of INTERMEDIATE products
# y***** = y(Tf*,Acou*,Aind*,a*,F*)
y11111 = Sn'*(x1 - sum(Z1,1))' + mf1111;
y10111 = Sn'*(eye(2340)-diag(sum(A011,1)))*L011*sum(F1,2) + mf1111;
y11001 = Sn'*(eye(2340)-diag(sum(A100,1)))*L100*sum(F1,2) + mf1111;
y10001 = Sn'*(eye(2340)-diag(sum(A0,1)))*L0*sum(F1,2) + mf1111;
y11101 = Sn'*(eye(2340)-diag(sum(A110,1)))*L110*sum(F1,2) + mf1111;
y10101 = Sn'*(eye(2340)-diag(sum(A010,1)))*L010*sum(F1,2) + mf1111;
y11011 = Sn'*(eye(2340)-diag(sum(A101,1)))*L101*sum(F1,2) + mf1111;
y10011 = Sn'*(eye(2340)-diag(sum(A001,1)))*L001*sum(F1,2) + mf1111;
y01110 = Sn'*(eye(2340)-diag(sum(A1,1)))*L1*sum(F0,2) + mf0000;
y00110 = Sn'*(eye(2340)-diag(sum(A011,1)))*L011*sum(F0,2) + mf0000;
y01000 = Sn'*(eye(2340)-diag(sum(A100,1)))*L100*sum(F0,2) + mf0000;
y00000 = Sn'*(x0 - sum(Z0,1))' + mf0000;
y01100 = Sn'*(eye(2340)-diag(sum(A110,1)))*L110*sum(F0,2) + mf0000;
y00100 = Sn'*(eye(2340)-diag(sum(A010,1)))*L010*sum(F0,2) + mf0000;
y01010 = Sn'*(eye(2340)-diag(sum(A101,1)))*L101*sum(F0,2) + mf0000;
y00010 = Sn'*(eye(2340)-diag(sum(A001,1)))*L001*sum(F0,2) + mf0000;
y11110 = Sn'*(eye(2340)-diag(sum(A1,1)))*L1*sum(F0,2) + mf1000;
y10110 = Sn'*(eye(2340)-diag(sum(A011,1)))*L011*sum(F0,2) + mf1000;
y11000 = Sn'*(eye(2340)-diag(sum(A100,1)))*L100*sum(F0,2) + mf1000;
y10000 = Sn'*(eye(2340)-diag(sum(A0,1)))*L0*sum(F0,2) + mf1000;
y11100 = Sn'*(eye(2340)-diag(sum(A110,1)))*L110*sum(F0,2) + mf1000;
y10100 = Sn'*(eye(2340)-diag(sum(A010,1)))*L010*sum(F0,2) + mf1000;
y11010 = Sn'*(eye(2340)-diag(sum(A101,1)))*L101*sum(F0,2) + mf1000;
y10010 = Sn'*(eye(2340)-diag(sum(A001,1)))*L001*sum(F0,2) + mf1000;
y01111 = Sn'*(eye(2340)-diag(sum(A1,1)))*L1*sum(F1,2) + mf0111;
y00111 = Sn'*(eye(2340)-diag(sum(A011,1)))*L011*sum(F1,2) + mf0111;
y01001 = Sn'*(eye(2340)-diag(sum(A100,1)))*L100*sum(F1,2) + mf0111;
y00001 = Sn'*(eye(2340)-diag(sum(A0,1)))*L0*sum(F1,2) + mf0111;
y01101 = Sn'*(eye(2340)-diag(sum(A110,1)))*L110*sum(F1,2) + mf0111;
y00101 = Sn'*(eye(2340)-diag(sum(A010,1)))*L010*sum(F1,2) + mf0111;
y01011 = Sn'*(eye(2340)-diag(sum(A101,1)))*L101*sum(F1,2) + mf0111;
y00011 = Sn'*(eye(2340)-diag(sum(A001,1)))*L001*sum(F1,2) + mf0111;

# clear unnecessary A*** and L*** terms to save space
clear A01* A10* L01* L10*
clear Z0 Z1

# print a message for monitoring the script execution (optional)
printf ("Now all data have been loaded, auxiliary matrices constructed.\n")

# run a bi-proportional balancing (equivalent to two RAS runs) to a matrix of ones – intermediate inputs
# r - partner (country of origin), s - importer (country in focus), i - industry of origin, j - producer (industry in focus)
# in the OECD ICIO 2018: Korea is No.19 and USA is No.36
for s = [19,36]
# prepare a matrix to store the results of the bi-proportinal balancing for each country s
    TCD_A = sparse(2340,2340);
    for j = 1:36
# for each s/jth column in dAcou_neg and dAcou_pos, generate a block-diagonal matrix with N K x K matrices of ones in its diagonal blocks; this implicitly involves transition from the country-industry to the industry-country classification, so each block in the matrix corresponds to product i
        TCD_A0 = kron(speye(36),ones(65));
# to apply the bi-proportional method of balancing, I adjust the block-diagonal matrix of ones TCD_A0 by scaling it row-wise and column-wise as in the first and second RAS runs (only two RAS runs are enough in this case!), so that it conforms with the row and column totals
# calculate a column vector of row multipliers
        rRAS = (P1*dAcou_pos(:,36*(s-1)+j))./sum(TCD_A0, 2);
# scale row-wise
        TCD_A1 = diag(rRAS)*TCD_A0;
# calculate a row vector of column multipliers
        sRAS = ((-P1*dAcou_neg(:,36*(s-1)+j))')./sum(TCD_A1, 1);
# address possible "NaN" in sRAS, replacing those with 0
		sRAS(isnan(sRAS)==1) = 0;
# scale column-wise
        TCD_A1 = TCD_A1*diag(sRAS);
# TCD_A1 is a block-diagonal matrix that shows the re-distribution between countries of origin (r) of each intermediate input (i) used in the production of the selected industry (j) in country (s)
# sequentially reorganize TCD_A1 into block columns and stack in TCD_A to finally obtain a KN x KN matrix for all industries j in country s
        for n = 1:36
            TCD_A(1+65*(n-1):65+65*(n-1),1+65*(j-1):65+65*(j-1)) = TCD_A1(1+65*(n-1):65+65*(n-1),1+65*(n-1):65+65*(n-1));
        endfor
# end loop for j
    endfor
    clear TCD_A1

# extract the information from TCD_A on trade creation (diversion, contraction) between country r and country s
# reshuffle rows and columns with the permutation matrix P, revert to the original country-industry classification
    TCD_A = P1'*TCD_A*P1;
    for r = [36,19]
# ignore a case where r = s
        if (r != s)
# construct a matrix of trade creation between countries s and r
			TC_A = sparse(2340,2340);
# copy the r-s block from TCD_A to the r-s block of TC_A
			TC_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1)) = TCD_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1));
# copy the r-s block of TCD_A to the s-s block of TC_A, changing sign
			TC_A(1+36*(s-1):36+36*(s-1),1+36*(s-1):36+36*(s-1)) = -TCD_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1));
# construct a matrix of trade diversion in between third countries and r
			TDin_A = sparse(2340,2340);
# block-transpose the row of r blocks from TCD_A, changing sign
			for t = 1:65
				TDin_A(1+36*(t-1):36+36*(t-1),1+36*(s-1):36+36*(s-1)) = -TCD_A(1+36*(r-1):36+36*(r-1),1+36*(t-1):36+36*(t-1));
			endfor
# calculate the TDin for country r as aggregate TDin in all r-t pairs less TC in the r-s pair
			TDin_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1)) = TCD_A(1+36*(r-1):36+36*(r-1),:)*Sk' - TCD_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1));
# set the s-s block to zero
			TDin_A(1+36*(s-1):36+36*(s-1),1+36*(s-1):36+36*(s-1)) = zeros(36);
# construct a matrix of trade diversion out between third countries and r
			TDout_A = sparse(2340,2340);
# copy the block column of r from TCD_A to the block column of s in TDout_A
			TDout_A(:,1+36*(s-1):36+36*(s-1)) = TCD_A(:,1+36*(r-1):36+36*(r-1));
# calculate TDout for country r as aggregate TDout in all t-r pairs less TR in the r-s pair, and change sign
			TDout_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1)) = -(Sk*TCD_A(:,1+36*(r-1):36+36*(r-1)) - TCD_A(1+36*(s-1):36+36*(s-1),1+36*(r-1):36+36*(r-1)));
# set the s-s block to zero
			TDout_A(1+36*(s-1):36+36*(s-1),1+36*(s-1):36+36*(s-1)) = zeros(36);
# construct a matrix of trade contraction (import substitution) between s and r
			TR_A = sparse(2340,2340);
# copy the s-r block from TCD_A to the s-s block of TR_A
			TR_A(1+36*(s-1):36+36*(s-1),1+36*(s-1):36+36*(s-1)) = TCD_A(1+36*(s-1):36+36*(s-1),1+36*(r-1):36+36*(r-1));
# copy the s-r block from TCD_A to the r-s block of TR_A, changing sign
			TR_A(1+36*(r-1):36+36*(r-1),1+36*(s-1):36+36*(s-1)) = -TCD_A(1+36*(s-1):36+36*(s-1),1+36*(r-1):36+36*(r-1));
                
# the rest of the (bilateral) component matrices for SDA - TC
			A1TCA11 = (Acou1-TC_A).*(Sk'*Aind1*diag(a1));
			A0TCA11 = (Acou0+TC_A).*(Sk'*Aind1*diag(a1));
			A1TCA00 = (Acou1-TC_A).*(Sk'*Aind0*diag(a0));
			A0TCA00 = (Acou0+TC_A).*(Sk'*Aind0*diag(a0));
			A1TCA10 = (Acou1-TC_A).*(Sk'*Aind1*diag(a0));
			A0TCA10 = (Acou0+TC_A).*(Sk'*Aind1*diag(a0));
			A1TCA01 = (Acou1-TC_A).*(Sk'*Aind0*diag(a1));
			A0TCA01 = (Acou0+TC_A).*(Sk'*Aind0*diag(a1));
			L1TCA11 = inv(eye(2340)-A1TCA11);
			L0TCA11 = inv(eye(2340)-A0TCA11);
			L1TCA00 = inv(eye(2340)-A1TCA00);
			L0TCA00 = inv(eye(2340)-A0TCA00);
			L1TCA10 = inv(eye(2340)-A1TCA10);
			L0TCA10 = inv(eye(2340)-A0TCA10);
			L1TCA01 = inv(eye(2340)-A1TCA01);
			L0TCA01 = inv(eye(2340)-A0TCA01);
# the rest 32 VALUE ADDED terms for changes in the country of orgin of INTERMEDIATE demand
# v**TCA*** = v(Tz*,Acou*,Aind*,a*,F*)
            v11TCA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA11),1)))*L1TCA11*sum(F1,2);
            v10TCA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA11),1)))*L0TCA11*sum(F1,2);
            v11TCA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA00),1)))*L1TCA00*sum(F1,2);
            v10TCA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA00),1)))*L0TCA00*sum(F1,2);
            v11TCA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA10),1)))*L1TCA10*sum(F1,2);
            v10TCA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA10),1)))*L0TCA10*sum(F1,2);
            v11TCA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA01),1)))*L1TCA01*sum(F1,2);
            v10TCA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA01),1)))*L0TCA01*sum(F1,2);
            v01TCA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA11),1)))*L1TCA11*sum(F0,2);
            v00TCA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA11),1)))*L0TCA11*sum(F0,2);
            v01TCA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA00),1)))*L1TCA00*sum(F0,2);
            v00TCA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA00),1)))*L0TCA00*sum(F0,2);
            v01TCA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA10),1)))*L1TCA10*sum(F0,2);
            v00TCA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA10),1)))*L0TCA10*sum(F0,2);
            v01TCA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA01),1)))*L1TCA01*sum(F0,2);
            v00TCA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA01),1)))*L0TCA01*sum(F0,2);
            v11TCA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA11),1)))*L1TCA11*sum(F0,2);
            v10TCA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA11),1)))*L0TCA11*sum(F0,2);
            v11TCA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA00),1)))*L1TCA00*sum(F0,2);
            v10TCA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA00),1)))*L0TCA00*sum(F0,2);
            v11TCA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA10),1)))*L1TCA10*sum(F0,2);
            v10TCA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA10),1)))*L0TCA10*sum(F0,2);
            v11TCA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TCA01),1)))*L1TCA01*sum(F0,2);
            v10TCA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TCA01),1)))*L0TCA01*sum(F0,2);
            v01TCA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA11),1)))*L1TCA11*sum(F1,2);
            v00TCA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA11),1)))*L0TCA11*sum(F1,2);
            v01TCA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA00),1)))*L1TCA00*sum(F1,2);
            v00TCA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA00),1)))*L0TCA00*sum(F1,2);
            v01TCA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA10),1)))*L1TCA10*sum(F1,2);
            v00TCA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA10),1)))*L0TCA10*sum(F1,2);
            v01TCA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TCA01),1)))*L1TCA01*sum(F1,2);
            v00TCA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TCA01),1)))*L0TCA01*sum(F1,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of INTERMEDIATE demand
# e**TCA*** = e(ec*,Acou*,Aind*,a*,F*)
            e11TCA111 = diag(ec1)*L1TCA11*sum(F1,2);
            e10TCA111 = diag(ec1)*L0TCA11*sum(F1,2);
            e11TCA001 = diag(ec1)*L1TCA00*sum(F1,2);
            e10TCA001 = diag(ec1)*L0TCA00*sum(F1,2);
            e11TCA101 = diag(ec1)*L1TCA10*sum(F1,2);
            e10TCA101 = diag(ec1)*L0TCA10*sum(F1,2);
            e11TCA011 = diag(ec1)*L1TCA01*sum(F1,2);
            e10TCA011 = diag(ec1)*L0TCA01*sum(F1,2);
            e01TCA110 = diag(ec0)*L1TCA11*sum(F0,2);
            e00TCA110 = diag(ec0)*L0TCA11*sum(F0,2);
            e01TCA000 = diag(ec0)*L1TCA00*sum(F0,2);
            e00TCA000 = diag(ec0)*L0TCA00*sum(F0,2);
            e01TCA100 = diag(ec0)*L1TCA10*sum(F0,2);
            e00TCA100 = diag(ec0)*L0TCA10*sum(F0,2);
            e01TCA010 = diag(ec0)*L1TCA01*sum(F0,2);
            e00TCA010 = diag(ec0)*L0TCA01*sum(F0,2);
            e11TCA110 = diag(ec1)*L1TCA11*sum(F0,2);
            e10TCA110 = diag(ec1)*L0TCA11*sum(F0,2);
            e11TCA000 = diag(ec1)*L1TCA00*sum(F0,2);
            e10TCA000 = diag(ec1)*L0TCA00*sum(F0,2);
            e11TCA100 = diag(ec1)*L1TCA10*sum(F0,2);
            e10TCA100 = diag(ec1)*L0TCA10*sum(F0,2);
            e11TCA010 = diag(ec1)*L1TCA01*sum(F0,2);
            e10TCA010 = diag(ec1)*L0TCA01*sum(F0,2);
            e01TCA111 = diag(ec0)*L1TCA11*sum(F1,2);
            e00TCA111 = diag(ec0)*L0TCA11*sum(F1,2);
            e01TCA001 = diag(ec0)*L1TCA00*sum(F1,2);
            e00TCA001 = diag(ec0)*L0TCA00*sum(F1,2);
            e01TCA101 = diag(ec0)*L1TCA10*sum(F1,2);
            e00TCA101 = diag(ec0)*L0TCA10*sum(F1,2);
            e01TCA011 = diag(ec0)*L1TCA01*sum(F1,2);
            e00TCA011 = diag(ec0)*L0TCA01*sum(F1,2);
# the rest 32 GDP terms for changes in the country of orgin of INTERMEDIATE demand
# y**TCA*** = y(Tf*,Acou*,Aind*,a*,F*)
            y11TCA111 = Sn'*(eye(2340)-diag(sum(A1TCA11,1)))*L1TCA11*sum(F1,2) + mf1111;
            y10TCA111 = Sn'*(eye(2340)-diag(sum(A0TCA11,1)))*L0TCA11*sum(F1,2) + mf1111;
            y11TCA001 = Sn'*(eye(2340)-diag(sum(A1TCA00,1)))*L1TCA00*sum(F1,2) + mf1111;
            y10TCA001 = Sn'*(eye(2340)-diag(sum(A0TCA00,1)))*L0TCA00*sum(F1,2) + mf1111;
            y11TCA101 = Sn'*(eye(2340)-diag(sum(A1TCA10,1)))*L1TCA10*sum(F1,2) + mf1111;
            y10TCA101 = Sn'*(eye(2340)-diag(sum(A0TCA10,1)))*L0TCA10*sum(F1,2) + mf1111;
            y11TCA011 = Sn'*(eye(2340)-diag(sum(A1TCA01,1)))*L1TCA01*sum(F1,2) + mf1111;
            y10TCA011 = Sn'*(eye(2340)-diag(sum(A0TCA01,1)))*L0TCA01*sum(F1,2) + mf1111;
            y01TCA110 = Sn'*(eye(2340)-diag(sum(A1TCA11,1)))*L1TCA11*sum(F0,2) + mf0000;
            y00TCA110 = Sn'*(eye(2340)-diag(sum(A0TCA11,1)))*L0TCA11*sum(F0,2) + mf0000;
            y01TCA000 = Sn'*(eye(2340)-diag(sum(A1TCA00,1)))*L1TCA00*sum(F0,2) + mf0000;
            y00TCA000 = Sn'*(eye(2340)-diag(sum(A0TCA00,1)))*L0TCA00*sum(F0,2) + mf0000;
            y01TCA100 = Sn'*(eye(2340)-diag(sum(A1TCA10,1)))*L1TCA10*sum(F0,2) + mf0000;
            y00TCA100 = Sn'*(eye(2340)-diag(sum(A0TCA10,1)))*L0TCA10*sum(F0,2) + mf0000;
            y01TCA010 = Sn'*(eye(2340)-diag(sum(A1TCA01,1)))*L1TCA01*sum(F0,2) + mf0000;
            y00TCA010 = Sn'*(eye(2340)-diag(sum(A0TCA01,1)))*L0TCA01*sum(F0,2) + mf0000;
            y11TCA110 = Sn'*(eye(2340)-diag(sum(A1TCA11,1)))*L1TCA11*sum(F0,2) + mf1000;
            y10TCA110 = Sn'*(eye(2340)-diag(sum(A0TCA11,1)))*L0TCA11*sum(F0,2) + mf1000;
            y11TCA000 = Sn'*(eye(2340)-diag(sum(A1TCA00,1)))*L1TCA00*sum(F0,2) + mf1000;
            y10TCA000 = Sn'*(eye(2340)-diag(sum(A0TCA00,1)))*L0TCA00*sum(F0,2) + mf1000;
            y11TCA100 = Sn'*(eye(2340)-diag(sum(A1TCA10,1)))*L1TCA10*sum(F0,2) + mf1000;
            y10TCA100 = Sn'*(eye(2340)-diag(sum(A0TCA10,1)))*L0TCA10*sum(F0,2) + mf1000;
            y11TCA010 = Sn'*(eye(2340)-diag(sum(A1TCA01,1)))*L1TCA01*sum(F0,2) + mf1000;
            y10TCA010 = Sn'*(eye(2340)-diag(sum(A0TCA01,1)))*L0TCA01*sum(F0,2) + mf1000;
            y01TCA111 = Sn'*(eye(2340)-diag(sum(A1TCA11,1)))*L1TCA11*sum(F1,2) + mf0111;
            y00TCA111 = Sn'*(eye(2340)-diag(sum(A0TCA11,1)))*L0TCA11*sum(F1,2) + mf0111;
            y01TCA001 = Sn'*(eye(2340)-diag(sum(A1TCA00,1)))*L1TCA00*sum(F1,2) + mf0111;
            y00TCA001 = Sn'*(eye(2340)-diag(sum(A0TCA00,1)))*L0TCA00*sum(F1,2) + mf0111;
            y01TCA101 = Sn'*(eye(2340)-diag(sum(A1TCA10,1)))*L1TCA10*sum(F1,2) + mf0111;
            y00TCA101 = Sn'*(eye(2340)-diag(sum(A0TCA10,1)))*L0TCA10*sum(F1,2) + mf0111;
            y01TCA011 = Sn'*(eye(2340)-diag(sum(A1TCA01,1)))*L1TCA01*sum(F1,2) + mf0111;
            y00TCA011 = Sn'*(eye(2340)-diag(sum(A0TCA01,1)))*L0TCA01*sum(F1,2) + mf0111;
# clear matrices that are not required until the next loop of r
            clear TC_A A0TCA* A1TCA* L0TCA* L1TCA*
# calculate the changes in VALUE ADDED as a weighted average of all TC-related terms - ADDITIVE
            dv_TC_A = (v1'-v11TCA111)/18 + (v10TCA111-v10111)/18 + (v11001-v11TCA001)/18 + (v10TCA001-v10001)/18 + (v11101-v11TCA101)/36 + (v10TCA101-v10101)/36 + (v11011-v11TCA011)/36 + (v10TCA011-v10011)/36 + (v01110-v01TCA110)/18 + (v00TCA110-v00110)/18 + (v01000-v01TCA000)/18 + (v00TCA000-v0')/18 + (v01100-v01TCA100)/36 + (v00TCA100-v00100)/36 + (v01010-v01TCA010)/36 + (v00TCA010-v00010)/36 + (v11110-v11TCA110)/36 + (v10TCA110-v10110)/36 + (v11000-v11TCA000)/36 + (v10TCA000-v10000)/36 + (v11100-v11TCA100)/72 + (v10TCA100-v10100)/72 + (v11010-v11TCA010)/72 + (v10TCA010-v10010)/72 + (v01111-v01TCA111)/36 + (v00TCA111-v00111)/36 + (v01001-v01TCA001)/36 + (v00TCA001-v00001)/36 + (v01101-v01TCA101)/72 + (v00TCA101-v00101)/72 + (v01011-v01TCA011)/72 + (v00TCA011-v00011)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TC-related terms - MULTIPLICATIVE
            rv_TC_A = (((v1'./v11TCA111) .* (v10TCA111./v10111) .* (v11001./v11TCA001) .* (v10TCA001./v10001)).^(1/18)) .* (((v11101./v11TCA101) .* (v10TCA101./v10101) .* (v11011./v11TCA011) .* (v10TCA011./v10011)).^(1/36)) .* (((v01110./v01TCA110) .* (v00TCA110./v00110) .* (v01000./v01TCA000) .* (v00TCA000./v0')).^(1/18)) .* (((v01100./v01TCA100) .* (v00TCA100./v00100) .* (v01010./v01TCA010) .* (v00TCA010./v00010)).^(1/36)) .* (((v11110./v11TCA110) .* (v10TCA110./v10110) .* (v11000./v11TCA000) .* (v10TCA000./v10000)).^(1/36)) .* (((v11100./v11TCA100) .* (v10TCA100./v10100) .* (v11010./v11TCA010) .* (v10TCA010./v10010)).^(1/72)) .* (((v01111./v01TCA111) .* (v00TCA111./v00111) .* (v01001./v01TCA001) .* (v00TCA001./v00001)).^(1/36)) .* (((v01101./v01TCA101) .* (v00TCA101./v00101) .* (v01011./v01TCA011) .* (v00TCA011./v00011)).^(1/72));
# NOTE: apparently, rv_TC_A and all rv_** may contain a few "NaN" entries. These, however, are unlikely to affect the multiplication of rv_**_A and rv_**_A later. "NaN" will only appear in those industries where zero is divided by zero.
# store results in structures
            dv.TC.A.(["r",num2str(r),"_s",num2str(s)]) = dv_TC_A;
            rv.TC.A.(["r",num2str(r),"_s",num2str(s)]) = rv_TC_A;
# calculate the changes in EMPLOYMENT as a weighted average of all TC-related terms - ADDITIVE
            de_TC_A = (e1-e11TCA111)/18 + (e10TCA111-e10111)/18 + (e11001-e11TCA001)/18 + (e10TCA001-e10001)/18 + (e11101-e11TCA101)/36 + (e10TCA101-e10101)/36 + (e11011-e11TCA011)/36 + (e10TCA011-e10011)/36 + (e01110-e01TCA110)/18 + (e00TCA110-e00110)/18 + (e01000-e01TCA000)/18 + (e00TCA000-e0)/18 + (e01100-e01TCA100)/36 + (e00TCA100-e00100)/36 + (e01010-e01TCA010)/36 + (e00TCA010-e00010)/36 + (e11110-e11TCA110)/36 + (e10TCA110-e10110)/36 + (e11000-e11TCA000)/36 + (e10TCA000-e10000)/36 + (e11100-e11TCA100)/72 + (e10TCA100-e10100)/72 + (e11010-e11TCA010)/72 + (e10TCA010-e10010)/72 + (e01111-e01TCA111)/36 + (e00TCA111-e00111)/36 + (e01001-e01TCA001)/36 + (e00TCA001-e00001)/36 + (e01101-e01TCA101)/72 + (e00TCA101-e00101)/72 + (e01011-e01TCA011)/72 + (e00TCA011-e00011)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TC-related terms - MULTIPLICATIVE
            re_TC_A = (((e1./e11TCA111) .* (e10TCA111./e10111) .* (e11001./e11TCA001) .* (e10TCA001./e10001)).^(1/18)) .* (((e11101./e11TCA101) .* (e10TCA101./e10101) .* (e11011./e11TCA011) .* (e10TCA011./e10011)).^(1/36)) .* (((e01110./e01TCA110) .* (e00TCA110./e00110) .* (e01000./e01TCA000) .* (e00TCA000./e0)).^(1/18)) .* (((e01100./e01TCA100) .* (e00TCA100./e00100) .* (e01010./e01TCA010) .* (e00TCA010./e00010)).^(1/36)) .* (((e11110./e11TCA110) .* (e10TCA110./e10110) .* (e11000./e11TCA000) .* (e10TCA000./e10000)).^(1/36)) .* (((e11100./e11TCA100) .* (e10TCA100./e10100) .* (e11010./e11TCA010) .* (e10TCA010./e10010)).^(1/72)) .* (((e01111./e01TCA111) .* (e00TCA111./e00111) .* (e01001./e01TCA001) .* (e00TCA001./e00001)).^(1/36)) .* (((e01101./e01TCA101) .* (e00TCA101./e00101) .* (e01011./e01TCA011) .* (e00TCA011./e00011)).^(1/72));
# store results in structures
            de.TC.A.(["r",num2str(r),"_s",num2str(s)]) = de_TC_A;
            re.TC.A.(["r",num2str(r),"_s",num2str(s)]) = re_TC_A;
# calculate the changes in GDP as a weighted average of all TC-related terms - ADDITIVE
            dy_TC_A = (y11111-y11TCA111)/18 + (y10TCA111-y10111)/18 + (y11001-y11TCA001)/18 + (y10TCA001-y10001)/18 + (y11101-y11TCA101)/36 + (y10TCA101-y10101)/36 + (y11011-y11TCA011)/36 + (y10TCA011-y10011)/36 + (y01110-y01TCA110)/18 + (y00TCA110-y00110)/18 + (y01000-y01TCA000)/18 + (y00TCA000-y00000)/18 + (y01100-y01TCA100)/36 + (y00TCA100-y00100)/36 + (y01010-y01TCA010)/36 + (y00TCA010-y00010)/36 + (y11110-y11TCA110)/36 + (y10TCA110-y10110)/36 + (y11000-y11TCA000)/36 + (y10TCA000-y10000)/36 + (y11100-y11TCA100)/72 + (y10TCA100-y10100)/72 + (y11010-y11TCA010)/72 + (y10TCA010-y10010)/72 + (y01111-y01TCA111)/36 + (y00TCA111-y00111)/36 + (y01001-y01TCA001)/36 + (y00TCA001-y00001)/36 + (y01101-y01TCA101)/72 + (y00TCA101-y00101)/72 + (y01011-y01TCA011)/72 + (y00TCA011-y00011)/72;
# calculate the changes in GDP as a weighted average of all TC-related terms - MULTIPLICATIVE
            ry_TC_A = (((y11111./y11TCA111) .* (y10TCA111./y10111) .* (y11001./y11TCA001) .* (y10TCA001./y10001)).^(1/18)) .* (((y11101./y11TCA101) .* (y10TCA101./y10101) .* (y11011./y11TCA011) .* (y10TCA011./y10011)).^(1/36)) .* (((y01110./y01TCA110) .* (y00TCA110./y00110) .* (y01000./y01TCA000) .* (y00TCA000./y00000)).^(1/18)) .* (((y01100./y01TCA100) .* (y00TCA100./y00100) .* (y01010./y01TCA010) .* (y00TCA010./y00010)).^(1/36)) .* (((y11110./y11TCA110) .* (y10TCA110./y10110) .* (y11000./y11TCA000) .* (y10TCA000./y10000)).^(1/36)) .* (((y11100./y11TCA100) .* (y10TCA100./y10100) .* (y11010./y11TCA010) .* (y10TCA010./y10010)).^(1/72)) .* (((y01111./y01TCA111) .* (y00TCA111./y00111) .* (y01001./y01TCA001) .* (y00TCA001./y00001)).^(1/36)) .* (((y01101./y01TCA101) .* (y00TCA101./y00101) .* (y01011./y01TCA011) .* (y00TCA011./y00011)).^(1/72));
# store results in structures
            dy.TC.A.(["r",num2str(r),"_s",num2str(s)]) = dy_TC_A;
            ry.TC.A.(["r",num2str(r),"_s",num2str(s)]) = ry_TC_A;

# the rest of the (bilateral) component matrices for SDA - TDin
			A1TDinA11 = (Acou1-TDin_A).*(Sk'*Aind1*diag(a1));
			A0TDinA11 = (Acou0+TDin_A).*(Sk'*Aind1*diag(a1));
			A1TDinA00 = (Acou1-TDin_A).*(Sk'*Aind0*diag(a0));
			A0TDinA00 = (Acou0+TDin_A).*(Sk'*Aind0*diag(a0));
			A1TDinA10 = (Acou1-TDin_A).*(Sk'*Aind1*diag(a0));
			A0TDinA10 = (Acou0+TDin_A).*(Sk'*Aind1*diag(a0));
			A1TDinA01 = (Acou1-TDin_A).*(Sk'*Aind0*diag(a1));
			A0TDinA01 = (Acou0+TDin_A).*(Sk'*Aind0*diag(a1));
			L1TDinA11 = inv(eye(2340)-A1TDinA11);
			L0TDinA11 = inv(eye(2340)-A0TDinA11);
			L1TDinA00 = inv(eye(2340)-A1TDinA00);
			L0TDinA00 = inv(eye(2340)-A0TDinA00);
			L1TDinA10 = inv(eye(2340)-A1TDinA10);
			L0TDinA10 = inv(eye(2340)-A0TDinA10);
			L1TDinA01 = inv(eye(2340)-A1TDinA01);
			L0TDinA01 = inv(eye(2340)-A0TDinA01);
# the rest 32 VALUE ADDED terms for changes in the country of orgin of INTERMEDIATE demand
# v**TDin*** = v(Tz*,Acou*,Aind*,a*,F*)
            v11TDinA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA11),1)))*L1TDinA11*sum(F1,2);
            v10TDinA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA11),1)))*L0TDinA11*sum(F1,2);
            v11TDinA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA00),1)))*L1TDinA00*sum(F1,2);
            v10TDinA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA00),1)))*L0TDinA00*sum(F1,2);
            v11TDinA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA10),1)))*L1TDinA10*sum(F1,2);
            v10TDinA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA10),1)))*L0TDinA10*sum(F1,2);
            v11TDinA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA01),1)))*L1TDinA01*sum(F1,2);
            v10TDinA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA01),1)))*L0TDinA01*sum(F1,2);
            v01TDinA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA11),1)))*L1TDinA11*sum(F0,2);
            v00TDinA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA11),1)))*L0TDinA11*sum(F0,2);
            v01TDinA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA00),1)))*L1TDinA00*sum(F0,2);
            v00TDinA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA00),1)))*L0TDinA00*sum(F0,2);
            v01TDinA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA10),1)))*L1TDinA10*sum(F0,2);
            v00TDinA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA10),1)))*L0TDinA10*sum(F0,2);
            v01TDinA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA01),1)))*L1TDinA01*sum(F0,2);
            v00TDinA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA01),1)))*L0TDinA01*sum(F0,2);
            v11TDinA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA11),1)))*L1TDinA11*sum(F0,2);
            v10TDinA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA11),1)))*L0TDinA11*sum(F0,2);
            v11TDinA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA00),1)))*L1TDinA00*sum(F0,2);
            v10TDinA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA00),1)))*L0TDinA00*sum(F0,2);
            v11TDinA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA10),1)))*L1TDinA10*sum(F0,2);
            v10TDinA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA10),1)))*L0TDinA10*sum(F0,2);
            v11TDinA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDinA01),1)))*L1TDinA01*sum(F0,2);
            v10TDinA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDinA01),1)))*L0TDinA01*sum(F0,2);
            v01TDinA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA11),1)))*L1TDinA11*sum(F1,2);
            v00TDinA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA11),1)))*L0TDinA11*sum(F1,2);
            v01TDinA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA00),1)))*L1TDinA00*sum(F1,2);
            v00TDinA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA00),1)))*L0TDinA00*sum(F1,2);
            v01TDinA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA10),1)))*L1TDinA10*sum(F1,2);
            v00TDinA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA10),1)))*L0TDinA10*sum(F1,2);
            v01TDinA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDinA01),1)))*L1TDinA01*sum(F1,2);
            v00TDinA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDinA01),1)))*L0TDinA01*sum(F1,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of INTERMEDIATE demand
# e**TDin*** = e(ec*,Acou*,Aind*,a*,F*)
            e11TDinA111 = diag(ec1)*L1TDinA11*sum(F1,2);
            e10TDinA111 = diag(ec1)*L0TDinA11*sum(F1,2);
            e11TDinA001 = diag(ec1)*L1TDinA00*sum(F1,2);
            e10TDinA001 = diag(ec1)*L0TDinA00*sum(F1,2);
            e11TDinA101 = diag(ec1)*L1TDinA10*sum(F1,2);
            e10TDinA101 = diag(ec1)*L0TDinA10*sum(F1,2);
            e11TDinA011 = diag(ec1)*L1TDinA01*sum(F1,2);
            e10TDinA011 = diag(ec1)*L0TDinA01*sum(F1,2);
            e01TDinA110 = diag(ec0)*L1TDinA11*sum(F0,2);
            e00TDinA110 = diag(ec0)*L0TDinA11*sum(F0,2);
            e01TDinA000 = diag(ec0)*L1TDinA00*sum(F0,2);
            e00TDinA000 = diag(ec0)*L0TDinA00*sum(F0,2);
            e01TDinA100 = diag(ec0)*L1TDinA10*sum(F0,2);
            e00TDinA100 = diag(ec0)*L0TDinA10*sum(F0,2);
            e01TDinA010 = diag(ec0)*L1TDinA01*sum(F0,2);
            e00TDinA010 = diag(ec0)*L0TDinA01*sum(F0,2);
            e11TDinA110 = diag(ec1)*L1TDinA11*sum(F0,2);
            e10TDinA110 = diag(ec1)*L0TDinA11*sum(F0,2);
            e11TDinA000 = diag(ec1)*L1TDinA00*sum(F0,2);
            e10TDinA000 = diag(ec1)*L0TDinA00*sum(F0,2);
            e11TDinA100 = diag(ec1)*L1TDinA10*sum(F0,2);
            e10TDinA100 = diag(ec1)*L0TDinA10*sum(F0,2);
            e11TDinA010 = diag(ec1)*L1TDinA01*sum(F0,2);
            e10TDinA010 = diag(ec1)*L0TDinA01*sum(F0,2);
            e01TDinA111 = diag(ec0)*L1TDinA11*sum(F1,2);
            e00TDinA111 = diag(ec0)*L0TDinA11*sum(F1,2);
            e01TDinA001 = diag(ec0)*L1TDinA00*sum(F1,2);
            e00TDinA001 = diag(ec0)*L0TDinA00*sum(F1,2);
            e01TDinA101 = diag(ec0)*L1TDinA10*sum(F1,2);
            e00TDinA101 = diag(ec0)*L0TDinA10*sum(F1,2);
            e01TDinA011 = diag(ec0)*L1TDinA01*sum(F1,2);
            e00TDinA011 = diag(ec0)*L0TDinA01*sum(F1,2);
# the rest 32 GDP terms for changes in the country of orgin of INTERMEDIATE demand
# y**TDin*** = y(Tf*,Acou*,Aind*,a*,F*)
            y11TDinA111 = Sn'*(eye(2340)-diag(sum(A1TDinA11,1)))*L1TDinA11*sum(F1,2) + mf1111;
            y10TDinA111 = Sn'*(eye(2340)-diag(sum(A0TDinA11,1)))*L0TDinA11*sum(F1,2) + mf1111;
            y11TDinA001 = Sn'*(eye(2340)-diag(sum(A1TDinA00,1)))*L1TDinA00*sum(F1,2) + mf1111;
            y10TDinA001 = Sn'*(eye(2340)-diag(sum(A0TDinA00,1)))*L0TDinA00*sum(F1,2) + mf1111;
            y11TDinA101 = Sn'*(eye(2340)-diag(sum(A1TDinA10,1)))*L1TDinA10*sum(F1,2) + mf1111;
            y10TDinA101 = Sn'*(eye(2340)-diag(sum(A0TDinA10,1)))*L0TDinA10*sum(F1,2) + mf1111;
            y11TDinA011 = Sn'*(eye(2340)-diag(sum(A1TDinA01,1)))*L1TDinA01*sum(F1,2) + mf1111;
            y10TDinA011 = Sn'*(eye(2340)-diag(sum(A0TDinA01,1)))*L0TDinA01*sum(F1,2) + mf1111;
            y01TDinA110 = Sn'*(eye(2340)-diag(sum(A1TDinA11,1)))*L1TDinA11*sum(F0,2) + mf0000;
            y00TDinA110 = Sn'*(eye(2340)-diag(sum(A0TDinA11,1)))*L0TDinA11*sum(F0,2) + mf0000;
            y01TDinA000 = Sn'*(eye(2340)-diag(sum(A1TDinA00,1)))*L1TDinA00*sum(F0,2) + mf0000;
            y00TDinA000 = Sn'*(eye(2340)-diag(sum(A0TDinA00,1)))*L0TDinA00*sum(F0,2) + mf0000;
            y01TDinA100 = Sn'*(eye(2340)-diag(sum(A1TDinA10,1)))*L1TDinA10*sum(F0,2) + mf0000;
            y00TDinA100 = Sn'*(eye(2340)-diag(sum(A0TDinA10,1)))*L0TDinA10*sum(F0,2) + mf0000;
            y01TDinA010 = Sn'*(eye(2340)-diag(sum(A1TDinA01,1)))*L1TDinA01*sum(F0,2) + mf0000;
            y00TDinA010 = Sn'*(eye(2340)-diag(sum(A0TDinA01,1)))*L0TDinA01*sum(F0,2) + mf0000;
            y11TDinA110 = Sn'*(eye(2340)-diag(sum(A1TDinA11,1)))*L1TDinA11*sum(F0,2) + mf1000;
            y10TDinA110 = Sn'*(eye(2340)-diag(sum(A0TDinA11,1)))*L0TDinA11*sum(F0,2) + mf1000;
            y11TDinA000 = Sn'*(eye(2340)-diag(sum(A1TDinA00,1)))*L1TDinA00*sum(F0,2) + mf1000;
            y10TDinA000 = Sn'*(eye(2340)-diag(sum(A0TDinA00,1)))*L0TDinA00*sum(F0,2) + mf1000;
            y11TDinA100 = Sn'*(eye(2340)-diag(sum(A1TDinA10,1)))*L1TDinA10*sum(F0,2) + mf1000;
            y10TDinA100 = Sn'*(eye(2340)-diag(sum(A0TDinA10,1)))*L0TDinA10*sum(F0,2) + mf1000;
            y11TDinA010 = Sn'*(eye(2340)-diag(sum(A1TDinA01,1)))*L1TDinA01*sum(F0,2) + mf1000;
            y10TDinA010 = Sn'*(eye(2340)-diag(sum(A0TDinA01,1)))*L0TDinA01*sum(F0,2) + mf1000;
            y01TDinA111 = Sn'*(eye(2340)-diag(sum(A1TDinA11,1)))*L1TDinA11*sum(F1,2) + mf0111;
            y00TDinA111 = Sn'*(eye(2340)-diag(sum(A0TDinA11,1)))*L0TDinA11*sum(F1,2) + mf0111;
            y01TDinA001 = Sn'*(eye(2340)-diag(sum(A1TDinA00,1)))*L1TDinA00*sum(F1,2) + mf0111;
            y00TDinA001 = Sn'*(eye(2340)-diag(sum(A0TDinA00,1)))*L0TDinA00*sum(F1,2) + mf0111;
            y01TDinA101 = Sn'*(eye(2340)-diag(sum(A1TDinA10,1)))*L1TDinA10*sum(F1,2) + mf0111;
            y00TDinA101 = Sn'*(eye(2340)-diag(sum(A0TDinA10,1)))*L0TDinA10*sum(F1,2) + mf0111;
            y01TDinA011 = Sn'*(eye(2340)-diag(sum(A1TDinA01,1)))*L1TDinA01*sum(F1,2) + mf0111;
            y00TDinA011 = Sn'*(eye(2340)-diag(sum(A0TDinA01,1)))*L0TDinA01*sum(F1,2) + mf0111;
# clear matrices that are not required until the next loop of r
            clear TDin_A A0TDinA* A1TDinA* L0TDinA* L1TDinA*
# calculate the changes in VALUE ADDED as a weighted average of all TDin-related terms - ADDITIVE
            dv_TDin_A = (v1'-v11TDinA111)/18 + (v10TDinA111-v10111)/18 + (v11001-v11TDinA001)/18 + (v10TDinA001-v10001)/18 + (v11101-v11TDinA101)/36 + (v10TDinA101-v10101)/36 + (v11011-v11TDinA011)/36 + (v10TDinA011-v10011)/36 + (v01110-v01TDinA110)/18 + (v00TDinA110-v00110)/18 + (v01000-v01TDinA000)/18 + (v00TDinA000-v0')/18 + (v01100-v01TDinA100)/36 + (v00TDinA100-v00100)/36 + (v01010-v01TDinA010)/36 + (v00TDinA010-v00010)/36 + (v11110-v11TDinA110)/36 + (v10TDinA110-v10110)/36 + (v11000-v11TDinA000)/36 + (v10TDinA000-v10000)/36 + (v11100-v11TDinA100)/72 + (v10TDinA100-v10100)/72 + (v11010-v11TDinA010)/72 + (v10TDinA010-v10010)/72 + (v01111-v01TDinA111)/36 + (v00TDinA111-v00111)/36 + (v01001-v01TDinA001)/36 + (v00TDinA001-v00001)/36 + (v01101-v01TDinA101)/72 + (v00TDinA101-v00101)/72 + (v01011-v01TDinA011)/72 + (v00TDinA011-v00011)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TDin-related terms - MULTIPLICATIVE
            rv_TDin_A = (((v1'./v11TDinA111) .* (v10TDinA111./v10111) .* (v11001./v11TDinA001) .* (v10TDinA001./v10001)).^(1/18)) .* (((v11101./v11TDinA101) .* (v10TDinA101./v10101) .* (v11011./v11TDinA011) .* (v10TDinA011./v10011)).^(1/36)) .* (((v01110./v01TDinA110) .* (v00TDinA110./v00110) .* (v01000./v01TDinA000) .* (v00TDinA000./v0')).^(1/18)) .* (((v01100./v01TDinA100) .* (v00TDinA100./v00100) .* (v01010./v01TDinA010) .* (v00TDinA010./v00010)).^(1/36)) .* (((v11110./v11TDinA110) .* (v10TDinA110./v10110) .* (v11000./v11TDinA000) .* (v10TDinA000./v10000)).^(1/36)) .* (((v11100./v11TDinA100) .* (v10TDinA100./v10100) .* (v11010./v11TDinA010) .* (v10TDinA010./v10010)).^(1/72)) .* (((v01111./v01TDinA111) .* (v00TDinA111./v00111) .* (v01001./v01TDinA001) .* (v00TDinA001./v00001)).^(1/36)) .* (((v01101./v01TDinA101) .* (v00TDinA101./v00101) .* (v01011./v01TDinA011) .* (v00TDinA011./v00011)).^(1/72));
# store results in structures
            dv.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = dv_TDin_A;
            rv.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = rv_TDin_A;
# calculate the changes in EMPLOYMENT as a weighted average of all TDin-related terms - ADDITIVE
            de_TDin_A = (e1-e11TDinA111)/18 + (e10TDinA111-e10111)/18 + (e11001-e11TDinA001)/18 + (e10TDinA001-e10001)/18 + (e11101-e11TDinA101)/36 + (e10TDinA101-e10101)/36 + (e11011-e11TDinA011)/36 + (e10TDinA011-e10011)/36 + (e01110-e01TDinA110)/18 + (e00TDinA110-e00110)/18 + (e01000-e01TDinA000)/18 + (e00TDinA000-e0)/18 + (e01100-e01TDinA100)/36 + (e00TDinA100-e00100)/36 + (e01010-e01TDinA010)/36 + (e00TDinA010-e00010)/36 + (e11110-e11TDinA110)/36 + (e10TDinA110-e10110)/36 + (e11000-e11TDinA000)/36 + (e10TDinA000-e10000)/36 + (e11100-e11TDinA100)/72 + (e10TDinA100-e10100)/72 + (e11010-e11TDinA010)/72 + (e10TDinA010-e10010)/72 + (e01111-e01TDinA111)/36 + (e00TDinA111-e00111)/36 + (e01001-e01TDinA001)/36 + (e00TDinA001-e00001)/36 + (e01101-e01TDinA101)/72 + (e00TDinA101-e00101)/72 + (e01011-e01TDinA011)/72 + (e00TDinA011-e00011)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TDin-related terms - MULTIPLICATIVE
            re_TDin_A = (((e1./e11TDinA111) .* (e10TDinA111./e10111) .* (e11001./e11TDinA001) .* (e10TDinA001./e10001)).^(1/18)) .* (((e11101./e11TDinA101) .* (e10TDinA101./e10101) .* (e11011./e11TDinA011) .* (e10TDinA011./e10011)).^(1/36)) .* (((e01110./e01TDinA110) .* (e00TDinA110./e00110) .* (e01000./e01TDinA000) .* (e00TDinA000./e0)).^(1/18)) .* (((e01100./e01TDinA100) .* (e00TDinA100./e00100) .* (e01010./e01TDinA010) .* (e00TDinA010./e00010)).^(1/36)) .* (((e11110./e11TDinA110) .* (e10TDinA110./e10110) .* (e11000./e11TDinA000) .* (e10TDinA000./e10000)).^(1/36)) .* (((e11100./e11TDinA100) .* (e10TDinA100./e10100) .* (e11010./e11TDinA010) .* (e10TDinA010./e10010)).^(1/72)) .* (((e01111./e01TDinA111) .* (e00TDinA111./e00111) .* (e01001./e01TDinA001) .* (e00TDinA001./e00001)).^(1/36)) .* (((e01101./e01TDinA101) .* (e00TDinA101./e00101) .* (e01011./e01TDinA011) .* (e00TDinA011./e00011)).^(1/72));
# store results in structures
            de.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = de_TDin_A;
            re.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = re_TDin_A;
# calculate the changes in GDP as a weighted average of all TDin-related terms - ADDITIVE
            dy_TDin_A = (y11111-y11TDinA111)/18 + (y10TDinA111-y10111)/18 + (y11001-y11TDinA001)/18 + (y10TDinA001-y10001)/18 + (y11101-y11TDinA101)/36 + (y10TDinA101-y10101)/36 + (y11011-y11TDinA011)/36 + (y10TDinA011-y10011)/36 + (y01110-y01TDinA110)/18 + (y00TDinA110-y00110)/18 + (y01000-y01TDinA000)/18 + (y00TDinA000-y00000)/18 + (y01100-y01TDinA100)/36 + (y00TDinA100-y00100)/36 + (y01010-y01TDinA010)/36 + (y00TDinA010-y00010)/36 + (y11110-y11TDinA110)/36 + (y10TDinA110-y10110)/36 + (y11000-y11TDinA000)/36 + (y10TDinA000-y10000)/36 + (y11100-y11TDinA100)/72 + (y10TDinA100-y10100)/72 + (y11010-y11TDinA010)/72 + (y10TDinA010-y10010)/72 + (y01111-y01TDinA111)/36 + (y00TDinA111-y00111)/36 + (y01001-y01TDinA001)/36 + (y00TDinA001-y00001)/36 + (y01101-y01TDinA101)/72 + (y00TDinA101-y00101)/72 + (y01011-y01TDinA011)/72 + (y00TDinA011-y00011)/72;
# calculate the changes in GDP as a weighted average of all TDin-related terms - MULTIPLICATIVE
            ry_TDin_A = (((y11111./y11TDinA111) .* (y10TDinA111./y10111) .* (y11001./y11TDinA001) .* (y10TDinA001./y10001)).^(1/18)) .* (((y11101./y11TDinA101) .* (y10TDinA101./y10101) .* (y11011./y11TDinA011) .* (y10TDinA011./y10011)).^(1/36)) .* (((y01110./y01TDinA110) .* (y00TDinA110./y00110) .* (y01000./y01TDinA000) .* (y00TDinA000./y00000)).^(1/18)) .* (((y01100./y01TDinA100) .* (y00TDinA100./y00100) .* (y01010./y01TDinA010) .* (y00TDinA010./y00010)).^(1/36)) .* (((y11110./y11TDinA110) .* (y10TDinA110./y10110) .* (y11000./y11TDinA000) .* (y10TDinA000./y10000)).^(1/36)) .* (((y11100./y11TDinA100) .* (y10TDinA100./y10100) .* (y11010./y11TDinA010) .* (y10TDinA010./y10010)).^(1/72)) .* (((y01111./y01TDinA111) .* (y00TDinA111./y00111) .* (y01001./y01TDinA001) .* (y00TDinA001./y00001)).^(1/36)) .* (((y01101./y01TDinA101) .* (y00TDinA101./y00101) .* (y01011./y01TDinA011) .* (y00TDinA011./y00011)).^(1/72));
# store results in structures
            dy.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = dy_TDin_A;
            ry.TDin.A.(["r",num2str(r),"_s",num2str(s)]) = ry_TDin_A;

# the rest of the (bilateral) component matrices for SDA - TDout
            A1TDoutA11 = (Acou1-TDout_A).*(Sk'*Aind1*diag(a1));
            A0TDoutA11 = (Acou0+TDout_A).*(Sk'*Aind1*diag(a1));
            A1TDoutA00 = (Acou1-TDout_A).*(Sk'*Aind0*diag(a0));
            A0TDoutA00 = (Acou0+TDout_A).*(Sk'*Aind0*diag(a0));
            A1TDoutA10 = (Acou1-TDout_A).*(Sk'*Aind1*diag(a0));
            A0TDoutA10 = (Acou0+TDout_A).*(Sk'*Aind1*diag(a0));
            A1TDoutA01 = (Acou1-TDout_A).*(Sk'*Aind0*diag(a1));
            A0TDoutA01 = (Acou0+TDout_A).*(Sk'*Aind0*diag(a1));
            L1TDoutA11 = inv(eye(2340)-A1TDoutA11);
            L0TDoutA11 = inv(eye(2340)-A0TDoutA11);
            L1TDoutA00 = inv(eye(2340)-A1TDoutA00);
            L0TDoutA00 = inv(eye(2340)-A0TDoutA00);
            L1TDoutA10 = inv(eye(2340)-A1TDoutA10);
            L0TDoutA10 = inv(eye(2340)-A0TDoutA10);
            L1TDoutA01 = inv(eye(2340)-A1TDoutA01);
            L0TDoutA01 = inv(eye(2340)-A0TDoutA01);
# the rest 32 VALUE ADDED terms for changes in the country of orgin of INTERMEDIATE demand
# v**TDout*** = v(Tz*,Acou*,Aind*,a*,F*)
            v11TDoutA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA11),1)))*L1TDoutA11*sum(F1,2);
            v10TDoutA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA11),1)))*L0TDoutA11*sum(F1,2);
            v11TDoutA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA00),1)))*L1TDoutA00*sum(F1,2);
            v10TDoutA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA00),1)))*L0TDoutA00*sum(F1,2);
            v11TDoutA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA10),1)))*L1TDoutA10*sum(F1,2);
            v10TDoutA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA10),1)))*L0TDoutA10*sum(F1,2);
            v11TDoutA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA01),1)))*L1TDoutA01*sum(F1,2);
            v10TDoutA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA01),1)))*L0TDoutA01*sum(F1,2);
            v01TDoutA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA11),1)))*L1TDoutA11*sum(F0,2);
            v00TDoutA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA11),1)))*L0TDoutA11*sum(F0,2);
            v01TDoutA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA00),1)))*L1TDoutA00*sum(F0,2);
            v00TDoutA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA00),1)))*L0TDoutA00*sum(F0,2);
            v01TDoutA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA10),1)))*L1TDoutA10*sum(F0,2);
            v00TDoutA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA10),1)))*L0TDoutA10*sum(F0,2);
            v01TDoutA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA01),1)))*L1TDoutA01*sum(F0,2);
            v00TDoutA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA01),1)))*L0TDoutA01*sum(F0,2);
            v11TDoutA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA11),1)))*L1TDoutA11*sum(F0,2);
            v10TDoutA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA11),1)))*L0TDoutA11*sum(F0,2);
            v11TDoutA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA00),1)))*L1TDoutA00*sum(F0,2);
            v10TDoutA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA00),1)))*L0TDoutA00*sum(F0,2);
            v11TDoutA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA10),1)))*L1TDoutA10*sum(F0,2);
            v10TDoutA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA10),1)))*L0TDoutA10*sum(F0,2);
            v11TDoutA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TDoutA01),1)))*L1TDoutA01*sum(F0,2);
            v10TDoutA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TDoutA01),1)))*L0TDoutA01*sum(F0,2);
            v01TDoutA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA11),1)))*L1TDoutA11*sum(F1,2);
            v00TDoutA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA11),1)))*L0TDoutA11*sum(F1,2);
            v01TDoutA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA00),1)))*L1TDoutA00*sum(F1,2);
            v00TDoutA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA00),1)))*L0TDoutA00*sum(F1,2);
            v01TDoutA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA10),1)))*L1TDoutA10*sum(F1,2);
            v00TDoutA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA10),1)))*L0TDoutA10*sum(F1,2);
            v01TDoutA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TDoutA01),1)))*L1TDoutA01*sum(F1,2);
            v00TDoutA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TDoutA01),1)))*L0TDoutA01*sum(F1,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of INTERMEDIATE demand
# e**TDout*** = e(ec*,Acou*,Aind*,a*,F*)
            e11TDoutA111 = diag(ec1)*L1TDoutA11*sum(F1,2);
            e10TDoutA111 = diag(ec1)*L0TDoutA11*sum(F1,2);
            e11TDoutA001 = diag(ec1)*L1TDoutA00*sum(F1,2);
            e10TDoutA001 = diag(ec1)*L0TDoutA00*sum(F1,2);
            e11TDoutA101 = diag(ec1)*L1TDoutA10*sum(F1,2);
            e10TDoutA101 = diag(ec1)*L0TDoutA10*sum(F1,2);
            e11TDoutA011 = diag(ec1)*L1TDoutA01*sum(F1,2);
            e10TDoutA011 = diag(ec1)*L0TDoutA01*sum(F1,2);
            e01TDoutA110 = diag(ec0)*L1TDoutA11*sum(F0,2);
            e00TDoutA110 = diag(ec0)*L0TDoutA11*sum(F0,2);
            e01TDoutA000 = diag(ec0)*L1TDoutA00*sum(F0,2);
            e00TDoutA000 = diag(ec0)*L0TDoutA00*sum(F0,2);
            e01TDoutA100 = diag(ec0)*L1TDoutA10*sum(F0,2);
            e00TDoutA100 = diag(ec0)*L0TDoutA10*sum(F0,2);
            e01TDoutA010 = diag(ec0)*L1TDoutA01*sum(F0,2);
            e00TDoutA010 = diag(ec0)*L0TDoutA01*sum(F0,2);
            e11TDoutA110 = diag(ec1)*L1TDoutA11*sum(F0,2);
            e10TDoutA110 = diag(ec1)*L0TDoutA11*sum(F0,2);
            e11TDoutA000 = diag(ec1)*L1TDoutA00*sum(F0,2);
            e10TDoutA000 = diag(ec1)*L0TDoutA00*sum(F0,2);
            e11TDoutA100 = diag(ec1)*L1TDoutA10*sum(F0,2);
            e10TDoutA100 = diag(ec1)*L0TDoutA10*sum(F0,2);
            e11TDoutA010 = diag(ec1)*L1TDoutA01*sum(F0,2);
            e10TDoutA010 = diag(ec1)*L0TDoutA01*sum(F0,2);
            e01TDoutA111 = diag(ec0)*L1TDoutA11*sum(F1,2);
            e00TDoutA111 = diag(ec0)*L0TDoutA11*sum(F1,2);
            e01TDoutA001 = diag(ec0)*L1TDoutA00*sum(F1,2);
            e00TDoutA001 = diag(ec0)*L0TDoutA00*sum(F1,2);
            e01TDoutA101 = diag(ec0)*L1TDoutA10*sum(F1,2);
            e00TDoutA101 = diag(ec0)*L0TDoutA10*sum(F1,2);
            e01TDoutA011 = diag(ec0)*L1TDoutA01*sum(F1,2);
            e00TDoutA011 = diag(ec0)*L0TDoutA01*sum(F1,2);
# the rest 32 GDP terms for changes in the country of orgin of INTERMEDIATE demand
# y**TDout** = y(Tf*,Acou*,Aind*,a*,F*)
            y11TDoutA111 = Sn'*(eye(2340)-diag(sum(A1TDoutA11,1)))*L1TDoutA11*sum(F1,2) + mf1111;
            y10TDoutA111 = Sn'*(eye(2340)-diag(sum(A0TDoutA11,1)))*L0TDoutA11*sum(F1,2) + mf1111;
            y11TDoutA001 = Sn'*(eye(2340)-diag(sum(A1TDoutA00,1)))*L1TDoutA00*sum(F1,2) + mf1111;
            y10TDoutA001 = Sn'*(eye(2340)-diag(sum(A0TDoutA00,1)))*L0TDoutA00*sum(F1,2) + mf1111;
            y11TDoutA101 = Sn'*(eye(2340)-diag(sum(A1TDoutA10,1)))*L1TDoutA10*sum(F1,2) + mf1111;
            y10TDoutA101 = Sn'*(eye(2340)-diag(sum(A0TDoutA10,1)))*L0TDoutA10*sum(F1,2) + mf1111;
            y11TDoutA011 = Sn'*(eye(2340)-diag(sum(A1TDoutA01,1)))*L1TDoutA01*sum(F1,2) + mf1111;
            y10TDoutA011 = Sn'*(eye(2340)-diag(sum(A0TDoutA01,1)))*L0TDoutA01*sum(F1,2) + mf1111;
            y01TDoutA110 = Sn'*(eye(2340)-diag(sum(A1TDoutA11,1)))*L1TDoutA11*sum(F0,2) + mf0000;
            y00TDoutA110 = Sn'*(eye(2340)-diag(sum(A0TDoutA11,1)))*L0TDoutA11*sum(F0,2) + mf0000;
            y01TDoutA000 = Sn'*(eye(2340)-diag(sum(A1TDoutA00,1)))*L1TDoutA00*sum(F0,2) + mf0000;
            y00TDoutA000 = Sn'*(eye(2340)-diag(sum(A0TDoutA00,1)))*L0TDoutA00*sum(F0,2) + mf0000;
            y01TDoutA100 = Sn'*(eye(2340)-diag(sum(A1TDoutA10,1)))*L1TDoutA10*sum(F0,2) + mf0000;
            y00TDoutA100 = Sn'*(eye(2340)-diag(sum(A0TDoutA10,1)))*L0TDoutA10*sum(F0,2) + mf0000;
            y01TDoutA010 = Sn'*(eye(2340)-diag(sum(A1TDoutA01,1)))*L1TDoutA01*sum(F0,2) + mf0000;
            y00TDoutA010 = Sn'*(eye(2340)-diag(sum(A0TDoutA01,1)))*L0TDoutA01*sum(F0,2) + mf0000;
            y11TDoutA110 = Sn'*(eye(2340)-diag(sum(A1TDoutA11,1)))*L1TDoutA11*sum(F0,2) + mf1000;
            y10TDoutA110 = Sn'*(eye(2340)-diag(sum(A0TDoutA11,1)))*L0TDoutA11*sum(F0,2) + mf1000;
            y11TDoutA000 = Sn'*(eye(2340)-diag(sum(A1TDoutA00,1)))*L1TDoutA00*sum(F0,2) + mf1000;
            y10TDoutA000 = Sn'*(eye(2340)-diag(sum(A0TDoutA00,1)))*L0TDoutA00*sum(F0,2) + mf1000;
            y11TDoutA100 = Sn'*(eye(2340)-diag(sum(A1TDoutA10,1)))*L1TDoutA10*sum(F0,2) + mf1000;
            y10TDoutA100 = Sn'*(eye(2340)-diag(sum(A0TDoutA10,1)))*L0TDoutA10*sum(F0,2) + mf1000;
            y11TDoutA010 = Sn'*(eye(2340)-diag(sum(A1TDoutA01,1)))*L1TDoutA01*sum(F0,2) + mf1000;
            y10TDoutA010 = Sn'*(eye(2340)-diag(sum(A0TDoutA01,1)))*L0TDoutA01*sum(F0,2) + mf1000;
            y01TDoutA111 = Sn'*(eye(2340)-diag(sum(A1TDoutA11,1)))*L1TDoutA11*sum(F1,2) + mf0111;
            y00TDoutA111 = Sn'*(eye(2340)-diag(sum(A0TDoutA11,1)))*L0TDoutA11*sum(F1,2) + mf0111;
            y01TDoutA001 = Sn'*(eye(2340)-diag(sum(A1TDoutA00,1)))*L1TDoutA00*sum(F1,2) + mf0111;
            y00TDoutA001 = Sn'*(eye(2340)-diag(sum(A0TDoutA00,1)))*L0TDoutA00*sum(F1,2) + mf0111;
            y01TDoutA101 = Sn'*(eye(2340)-diag(sum(A1TDoutA10,1)))*L1TDoutA10*sum(F1,2) + mf0111;
            y00TDoutA101 = Sn'*(eye(2340)-diag(sum(A0TDoutA10,1)))*L0TDoutA10*sum(F1,2) + mf0111;
            y01TDoutA011 = Sn'*(eye(2340)-diag(sum(A1TDoutA01,1)))*L1TDoutA01*sum(F1,2) + mf0111;
            y00TDoutA011 = Sn'*(eye(2340)-diag(sum(A0TDoutA01,1)))*L0TDoutA01*sum(F1,2) + mf0111;
# clear matrices that are not required until the next loop of r
            clear TDout_A A0TDoutA* A1TDoutA* L0TDoutA* L1TDoutA*
# calculate the changes in VALUE ADDED as a weighted average of all TDout-related terms - ADDITIVE
            dv_TDout_A = (v1'-v11TDoutA111)/18 + (v10TDoutA111-v10111)/18 + (v11001-v11TDoutA001)/18 + (v10TDoutA001-v10001)/18 + (v11101-v11TDoutA101)/36 + (v10TDoutA101-v10101)/36 + (v11011-v11TDoutA011)/36 + (v10TDoutA011-v10011)/36 + (v01110-v01TDoutA110)/18 + (v00TDoutA110-v00110)/18 + (v01000-v01TDoutA000)/18 + (v00TDoutA000-v0')/18 + (v01100-v01TDoutA100)/36 + (v00TDoutA100-v00100)/36 + (v01010-v01TDoutA010)/36 + (v00TDoutA010-v00010)/36 + (v11110-v11TDoutA110)/36 + (v10TDoutA110-v10110)/36 + (v11000-v11TDoutA000)/36 + (v10TDoutA000-v10000)/36 + (v11100-v11TDoutA100)/72 + (v10TDoutA100-v10100)/72 + (v11010-v11TDoutA010)/72 + (v10TDoutA010-v10010)/72 + (v01111-v01TDoutA111)/36 + (v00TDoutA111-v00111)/36 + (v01001-v01TDoutA001)/36 + (v00TDoutA001-v00001)/36 + (v01101-v01TDoutA101)/72 + (v00TDoutA101-v00101)/72 + (v01011-v01TDoutA011)/72 + (v00TDoutA011-v00011)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TDout-related terms - MULTIPLICATIVE
            rv_TDout_A = (((v1'./v11TDoutA111) .* (v10TDoutA111./v10111) .* (v11001./v11TDoutA001) .* (v10TDoutA001./v10001)).^(1/18)) .* (((v11101./v11TDoutA101) .* (v10TDoutA101./v10101) .* (v11011./v11TDoutA011) .* (v10TDoutA011./v10011)).^(1/36)) .* (((v01110./v01TDoutA110) .* (v00TDoutA110./v00110) .* (v01000./v01TDoutA000) .* (v00TDoutA000./v0')).^(1/18)) .* (((v01100./v01TDoutA100) .* (v00TDoutA100./v00100) .* (v01010./v01TDoutA010) .* (v00TDoutA010./v00010)).^(1/36)) .* (((v11110./v11TDoutA110) .* (v10TDoutA110./v10110) .* (v11000./v11TDoutA000) .* (v10TDoutA000./v10000)).^(1/36)) .* (((v11100./v11TDoutA100) .* (v10TDoutA100./v10100) .* (v11010./v11TDoutA010) .* (v10TDoutA010./v10010)).^(1/72)) .* (((v01111./v01TDoutA111) .* (v00TDoutA111./v00111) .* (v01001./v01TDoutA001) .* (v00TDoutA001./v00001)).^(1/36)) .* (((v01101./v01TDoutA101) .* (v00TDoutA101./v00101) .* (v01011./v01TDoutA011) .* (v00TDoutA011./v00011)).^(1/72));
# store results in structures
            dv.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = dv_TDout_A;
            rv.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = rv_TDout_A;
# calculate the changes in EMPLOYMENT as a weighted average of all TDout-related terms - ADDITIVE
            de_TDout_A = (e1-e11TDoutA111)/18 + (e10TDoutA111-e10111)/18 + (e11001-e11TDoutA001)/18 + (e10TDoutA001-e10001)/18 + (e11101-e11TDoutA101)/36 + (e10TDoutA101-e10101)/36 + (e11011-e11TDoutA011)/36 + (e10TDoutA011-e10011)/36 + (e01110-e01TDoutA110)/18 + (e00TDoutA110-e00110)/18 + (e01000-e01TDoutA000)/18 + (e00TDoutA000-e0)/18 + (e01100-e01TDoutA100)/36 + (e00TDoutA100-e00100)/36 + (e01010-e01TDoutA010)/36 + (e00TDoutA010-e00010)/36 + (e11110-e11TDoutA110)/36 + (e10TDoutA110-e10110)/36 + (e11000-e11TDoutA000)/36 + (e10TDoutA000-e10000)/36 + (e11100-e11TDoutA100)/72 + (e10TDoutA100-e10100)/72 + (e11010-e11TDoutA010)/72 + (e10TDoutA010-e10010)/72 + (e01111-e01TDoutA111)/36 + (e00TDoutA111-e00111)/36 + (e01001-e01TDoutA001)/36 + (e00TDoutA001-e00001)/36 + (e01101-e01TDoutA101)/72 + (e00TDoutA101-e00101)/72 + (e01011-e01TDoutA011)/72 + (e00TDoutA011-e00011)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TDout-related terms - MULTIPLICATIVE
            re_TDout_A = (((e1./e11TDoutA111) .* (e10TDoutA111./e10111) .* (e11001./e11TDoutA001) .* (e10TDoutA001./e10001)).^(1/18)) .* (((e11101./e11TDoutA101) .* (e10TDoutA101./e10101) .* (e11011./e11TDoutA011) .* (e10TDoutA011./e10011)).^(1/36)) .* (((e01110./e01TDoutA110) .* (e00TDoutA110./e00110) .* (e01000./e01TDoutA000) .* (e00TDoutA000./e0)).^(1/18)) .* (((e01100./e01TDoutA100) .* (e00TDoutA100./e00100) .* (e01010./e01TDoutA010) .* (e00TDoutA010./e00010)).^(1/36)) .* (((e11110./e11TDoutA110) .* (e10TDoutA110./e10110) .* (e11000./e11TDoutA000) .* (e10TDoutA000./e10000)).^(1/36)) .* (((e11100./e11TDoutA100) .* (e10TDoutA100./e10100) .* (e11010./e11TDoutA010) .* (e10TDoutA010./e10010)).^(1/72)) .* (((e01111./e01TDoutA111) .* (e00TDoutA111./e00111) .* (e01001./e01TDoutA001) .* (e00TDoutA001./e00001)).^(1/36)) .* (((e01101./e01TDoutA101) .* (e00TDoutA101./e00101) .* (e01011./e01TDoutA011) .* (e00TDoutA011./e00011)).^(1/72));
# store results in structures
            de.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = de_TDout_A;
            re.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = re_TDout_A;
# calculate the changes in GDP as a weighted average of all TDout-related terms - ADDITIVE
            dy_TDout_A = (y11111-y11TDoutA111)/18 + (y10TDoutA111-y10111)/18 + (y11001-y11TDoutA001)/18 + (y10TDoutA001-y10001)/18 + (y11101-y11TDoutA101)/36 + (y10TDoutA101-y10101)/36 + (y11011-y11TDoutA011)/36 + (y10TDoutA011-y10011)/36 + (y01110-y01TDoutA110)/18 + (y00TDoutA110-y00110)/18 + (y01000-y01TDoutA000)/18 + (y00TDoutA000-y00000)/18 + (y01100-y01TDoutA100)/36 + (y00TDoutA100-y00100)/36 + (y01010-y01TDoutA010)/36 + (y00TDoutA010-y00010)/36 + (y11110-y11TDoutA110)/36 + (y10TDoutA110-y10110)/36 + (y11000-y11TDoutA000)/36 + (y10TDoutA000-y10000)/36 + (y11100-y11TDoutA100)/72 + (y10TDoutA100-y10100)/72 + (y11010-y11TDoutA010)/72 + (y10TDoutA010-y10010)/72 + (y01111-y01TDoutA111)/36 + (y00TDoutA111-y00111)/36 + (y01001-y01TDoutA001)/36 + (y00TDoutA001-y00001)/36 + (y01101-y01TDoutA101)/72 + (y00TDoutA101-y00101)/72 + (y01011-y01TDoutA011)/72 + (y00TDoutA011-y00011)/72;
# calculate the changes in GDP as a weighted average of all TDout-related terms - MULTIPLICATIVE
            ry_TDout_A = (((y11111./y11TDoutA111) .* (y10TDoutA111./y10111) .* (y11001./y11TDoutA001) .* (y10TDoutA001./y10001)).^(1/18)) .* (((y11101./y11TDoutA101) .* (y10TDoutA101./y10101) .* (y11011./y11TDoutA011) .* (y10TDoutA011./y10011)).^(1/36)) .* (((y01110./y01TDoutA110) .* (y00TDoutA110./y00110) .* (y01000./y01TDoutA000) .* (y00TDoutA000./y00000)).^(1/18)) .* (((y01100./y01TDoutA100) .* (y00TDoutA100./y00100) .* (y01010./y01TDoutA010) .* (y00TDoutA010./y00010)).^(1/36)) .* (((y11110./y11TDoutA110) .* (y10TDoutA110./y10110) .* (y11000./y11TDoutA000) .* (y10TDoutA000./y10000)).^(1/36)) .* (((y11100./y11TDoutA100) .* (y10TDoutA100./y10100) .* (y11010./y11TDoutA010) .* (y10TDoutA010./y10010)).^(1/72)) .* (((y01111./y01TDoutA111) .* (y00TDoutA111./y00111) .* (y01001./y01TDoutA001) .* (y00TDoutA001./y00001)).^(1/36)) .* (((y01101./y01TDoutA101) .* (y00TDoutA101./y00101) .* (y01011./y01TDoutA011) .* (y00TDoutA011./y00011)).^(1/72));
# store results in structures
            dy.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = dy_TDout_A;
            ry.TDout.A.(["r",num2str(r),"_s",num2str(s)]) = ry_TDout_A;

# the rest of the (bilateral) component matrices for SDA - TR
			A1TRA11 = (Acou1-TR_A).*(Sk'*Aind1*diag(a1));
			A0TRA11 = (Acou0+TR_A).*(Sk'*Aind1*diag(a1));
			A1TRA00 = (Acou1-TR_A).*(Sk'*Aind0*diag(a0));
			A0TRA00 = (Acou0+TR_A).*(Sk'*Aind0*diag(a0));
			A1TRA10 = (Acou1-TR_A).*(Sk'*Aind1*diag(a0));
			A0TRA10 = (Acou0+TR_A).*(Sk'*Aind1*diag(a0));
			A1TRA01 = (Acou1-TR_A).*(Sk'*Aind0*diag(a1));
			A0TRA01 = (Acou0+TR_A).*(Sk'*Aind0*diag(a1));
			L1TRA11 = inv(eye(2340)-A1TRA11);
			L0TRA11 = inv(eye(2340)-A0TRA11);
			L1TRA00 = inv(eye(2340)-A1TRA00);
			L0TRA00 = inv(eye(2340)-A0TRA00);
			L1TRA10 = inv(eye(2340)-A1TRA10);
			L0TRA10 = inv(eye(2340)-A0TRA10);
			L1TRA01 = inv(eye(2340)-A1TRA01);
			L0TRA01 = inv(eye(2340)-A0TRA01);
# the rest 32 VALUE ADDED terms for changes in the country of orgin of INTERMEDIATE demand
# v**TR*** = v(Tz*,Acou*,Aind*,a*,F*)
			v11TRA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA11),1)))*L1TRA11*sum(F1,2);
			v10TRA111 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA11),1)))*L0TRA11*sum(F1,2);
			v11TRA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA00),1)))*L1TRA00*sum(F1,2);
			v10TRA001 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA00),1)))*L0TRA00*sum(F1,2);
			v11TRA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA10),1)))*L1TRA10*sum(F1,2);
			v10TRA101 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA10),1)))*L0TRA10*sum(F1,2);
			v11TRA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA01),1)))*L1TRA01*sum(F1,2);
			v10TRA011 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA01),1)))*L0TRA01*sum(F1,2);
			v01TRA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA11),1)))*L1TRA11*sum(F0,2);
			v00TRA110 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA11),1)))*L0TRA11*sum(F0,2);
			v01TRA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA00),1)))*L1TRA00*sum(F0,2);
			v00TRA000 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA00),1)))*L0TRA00*sum(F0,2);
			v01TRA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA10),1)))*L1TRA10*sum(F0,2);
			v00TRA100 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA10),1)))*L0TRA10*sum(F0,2);
			v01TRA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA01),1)))*L1TRA01*sum(F0,2);
			v00TRA010 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA01),1)))*L0TRA01*sum(F0,2);
			v11TRA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA11),1)))*L1TRA11*sum(F0,2);
			v10TRA110 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA11),1)))*L0TRA11*sum(F0,2);
			v11TRA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA00),1)))*L1TRA00*sum(F0,2);
			v10TRA000 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA00),1)))*L0TRA00*sum(F0,2);
			v11TRA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA10),1)))*L1TRA10*sum(F0,2);
			v10TRA100 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA10),1)))*L0TRA10*sum(F0,2);
			v11TRA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1TRA01),1)))*L1TRA01*sum(F0,2);
			v10TRA010 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0TRA01),1)))*L0TRA01*sum(F0,2);
			v01TRA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA11),1)))*L1TRA11*sum(F1,2);
			v00TRA111 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA11),1)))*L0TRA11*sum(F1,2);
			v01TRA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA00),1)))*L1TRA00*sum(F1,2);
			v00TRA001 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA00),1)))*L0TRA00*sum(F1,2);
			v01TRA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA10),1)))*L1TRA10*sum(F1,2);
			v00TRA101 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA10),1)))*L0TRA10*sum(F1,2);
			v01TRA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1TRA01),1)))*L1TRA01*sum(F1,2);
			v00TRA011 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0TRA01),1)))*L0TRA01*sum(F1,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of INTERMEDIATE demand
# e**TR*** = e(ec*,Acou*,Aind*,a*,F*)
            e11TRA111 = diag(ec1)*L1TRA11*sum(F1,2);
            e10TRA111 = diag(ec1)*L0TRA11*sum(F1,2);
            e11TRA001 = diag(ec1)*L1TRA00*sum(F1,2);
            e10TRA001 = diag(ec1)*L0TRA00*sum(F1,2);
            e11TRA101 = diag(ec1)*L1TRA10*sum(F1,2);
            e10TRA101 = diag(ec1)*L0TRA10*sum(F1,2);
            e11TRA011 = diag(ec1)*L1TRA01*sum(F1,2);
            e10TRA011 = diag(ec1)*L0TRA01*sum(F1,2);
            e01TRA110 = diag(ec0)*L1TRA11*sum(F0,2);
            e00TRA110 = diag(ec0)*L0TRA11*sum(F0,2);
            e01TRA000 = diag(ec0)*L1TRA00*sum(F0,2);
            e00TRA000 = diag(ec0)*L0TRA00*sum(F0,2);
            e01TRA100 = diag(ec0)*L1TRA10*sum(F0,2);
            e00TRA100 = diag(ec0)*L0TRA10*sum(F0,2);
            e01TRA010 = diag(ec0)*L1TRA01*sum(F0,2);
            e00TRA010 = diag(ec0)*L0TRA01*sum(F0,2);
            e11TRA110 = diag(ec1)*L1TRA11*sum(F0,2);
            e10TRA110 = diag(ec1)*L0TRA11*sum(F0,2);
            e11TRA000 = diag(ec1)*L1TRA00*sum(F0,2);
            e10TRA000 = diag(ec1)*L0TRA00*sum(F0,2);
            e11TRA100 = diag(ec1)*L1TRA10*sum(F0,2);
            e10TRA100 = diag(ec1)*L0TRA10*sum(F0,2);
            e11TRA010 = diag(ec1)*L1TRA01*sum(F0,2);
            e10TRA010 = diag(ec1)*L0TRA01*sum(F0,2);
            e01TRA111 = diag(ec0)*L1TRA11*sum(F1,2);
            e00TRA111 = diag(ec0)*L0TRA11*sum(F1,2);
            e01TRA001 = diag(ec0)*L1TRA00*sum(F1,2);
            e00TRA001 = diag(ec0)*L0TRA00*sum(F1,2);
            e01TRA101 = diag(ec0)*L1TRA10*sum(F1,2);
            e00TRA101 = diag(ec0)*L0TRA10*sum(F1,2);
            e01TRA011 = diag(ec0)*L1TRA01*sum(F1,2);
            e00TRA011 = diag(ec0)*L0TRA01*sum(F1,2);
# the rest 32 GDP terms for changes in the country of orgin of INTERMEDIATE demand
# y**TR*** = y(Tf*,Acou*,Aind*,a*,F*)
            y11TRA111 = Sn'*(eye(2340)-diag(sum(A1TRA11,1)))*L1TRA11*sum(F1,2) + mf1111;
            y10TRA111 = Sn'*(eye(2340)-diag(sum(A0TRA11,1)))*L0TRA11*sum(F1,2) + mf1111;
            y11TRA001 = Sn'*(eye(2340)-diag(sum(A1TRA00,1)))*L1TRA00*sum(F1,2) + mf1111;
            y10TRA001 = Sn'*(eye(2340)-diag(sum(A0TRA00,1)))*L0TRA00*sum(F1,2) + mf1111;
            y11TRA101 = Sn'*(eye(2340)-diag(sum(A1TRA10,1)))*L1TRA10*sum(F1,2) + mf1111;
            y10TRA101 = Sn'*(eye(2340)-diag(sum(A0TRA10,1)))*L0TRA10*sum(F1,2) + mf1111;
            y11TRA011 = Sn'*(eye(2340)-diag(sum(A1TRA01,1)))*L1TRA01*sum(F1,2) + mf1111;
            y10TRA011 = Sn'*(eye(2340)-diag(sum(A0TRA01,1)))*L0TRA01*sum(F1,2) + mf1111;
            y01TRA110 = Sn'*(eye(2340)-diag(sum(A1TRA11,1)))*L1TRA11*sum(F0,2) + mf0000;
            y00TRA110 = Sn'*(eye(2340)-diag(sum(A0TRA11,1)))*L0TRA11*sum(F0,2) + mf0000;
            y01TRA000 = Sn'*(eye(2340)-diag(sum(A1TRA00,1)))*L1TRA00*sum(F0,2) + mf0000;
            y00TRA000 = Sn'*(eye(2340)-diag(sum(A0TRA00,1)))*L0TRA00*sum(F0,2) + mf0000;
            y01TRA100 = Sn'*(eye(2340)-diag(sum(A1TRA10,1)))*L1TRA10*sum(F0,2) + mf0000;
            y00TRA100 = Sn'*(eye(2340)-diag(sum(A0TRA10,1)))*L0TRA10*sum(F0,2) + mf0000;
            y01TRA010 = Sn'*(eye(2340)-diag(sum(A1TRA01,1)))*L1TRA01*sum(F0,2) + mf0000;
            y00TRA010 = Sn'*(eye(2340)-diag(sum(A0TRA01,1)))*L0TRA01*sum(F0,2) + mf0000;
            y11TRA110 = Sn'*(eye(2340)-diag(sum(A1TRA11,1)))*L1TRA11*sum(F0,2) + mf1000;
            y10TRA110 = Sn'*(eye(2340)-diag(sum(A0TRA11,1)))*L0TRA11*sum(F0,2) + mf1000;
            y11TRA000 = Sn'*(eye(2340)-diag(sum(A1TRA00,1)))*L1TRA00*sum(F0,2) + mf1000;
            y10TRA000 = Sn'*(eye(2340)-diag(sum(A0TRA00,1)))*L0TRA00*sum(F0,2) + mf1000;
            y11TRA100 = Sn'*(eye(2340)-diag(sum(A1TRA10,1)))*L1TRA10*sum(F0,2) + mf1000;
            y10TRA100 = Sn'*(eye(2340)-diag(sum(A0TRA10,1)))*L0TRA10*sum(F0,2) + mf1000;
            y11TRA010 = Sn'*(eye(2340)-diag(sum(A1TRA01,1)))*L1TRA01*sum(F0,2) + mf1000;
            y10TRA010 = Sn'*(eye(2340)-diag(sum(A0TRA01,1)))*L0TRA01*sum(F0,2) + mf1000;
            y01TRA111 = Sn'*(eye(2340)-diag(sum(A1TRA11,1)))*L1TRA11*sum(F1,2) + mf0111;
            y00TRA111 = Sn'*(eye(2340)-diag(sum(A0TRA11,1)))*L0TRA11*sum(F1,2) + mf0111;
            y01TRA001 = Sn'*(eye(2340)-diag(sum(A1TRA00,1)))*L1TRA00*sum(F1,2) + mf0111;
            y00TRA001 = Sn'*(eye(2340)-diag(sum(A0TRA00,1)))*L0TRA00*sum(F1,2) + mf0111;
            y01TRA101 = Sn'*(eye(2340)-diag(sum(A1TRA10,1)))*L1TRA10*sum(F1,2) + mf0111;
            y00TRA101 = Sn'*(eye(2340)-diag(sum(A0TRA10,1)))*L0TRA10*sum(F1,2) + mf0111;
            y01TRA011 = Sn'*(eye(2340)-diag(sum(A1TRA01,1)))*L1TRA01*sum(F1,2) + mf0111;
            y00TRA011 = Sn'*(eye(2340)-diag(sum(A0TRA01,1)))*L0TRA01*sum(F1,2) + mf0111;
# clear matrices that are not required until the next loop of r
            clear TR_A A0TRA* A1TRA* L0TRA* L1TRA*
# calculate the changes in VALUE ADDED as a weighted average of all TR-related terms - ADDITIVE
            dv_TR_A = (v1'-v11TRA111)/18 + (v10TRA111-v10111)/18 + (v11001-v11TRA001)/18 + (v10TRA001-v10001)/18 + (v11101-v11TRA101)/36 + (v10TRA101-v10101)/36 + (v11011-v11TRA011)/36 + (v10TRA011-v10011)/36 + (v01110-v01TRA110)/18 + (v00TRA110-v00110)/18 + (v01000-v01TRA000)/18 + (v00TRA000-v0')/18 + (v01100-v01TRA100)/36 + (v00TRA100-v00100)/36 + (v01010-v01TRA010)/36 + (v00TRA010-v00010)/36 + (v11110-v11TRA110)/36 + (v10TRA110-v10110)/36 + (v11000-v11TRA000)/36 + (v10TRA000-v10000)/36 + (v11100-v11TRA100)/72 + (v10TRA100-v10100)/72 + (v11010-v11TRA010)/72 + (v10TRA010-v10010)/72 + (v01111-v01TRA111)/36 + (v00TRA111-v00111)/36 + (v01001-v01TRA001)/36 + (v00TRA001-v00001)/36 + (v01101-v01TRA101)/72 + (v00TRA101-v00101)/72 + (v01011-v01TRA011)/72 + (v00TRA011-v00011)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TR-related terms - MULTIPLICATIVE
            rv_TR_A = (((v1'./v11TRA111) .* (v10TRA111./v10111) .* (v11001./v11TRA001) .* (v10TRA001./v10001)).^(1/18)) .* (((v11101./v11TRA101) .* (v10TRA101./v10101) .* (v11011./v11TRA011) .* (v10TRA011./v10011)).^(1/36)) .* (((v01110./v01TRA110) .* (v00TRA110./v00110) .* (v01000./v01TRA000) .* (v00TRA000./v0')).^(1/18)) .* (((v01100./v01TRA100) .* (v00TRA100./v00100) .* (v01010./v01TRA010) .* (v00TRA010./v00010)).^(1/36)) .* (((v11110./v11TRA110) .* (v10TRA110./v10110) .* (v11000./v11TRA000) .* (v10TRA000./v10000)).^(1/36)) .* (((v11100./v11TRA100) .* (v10TRA100./v10100) .* (v11010./v11TRA010) .* (v10TRA010./v10010)).^(1/72)) .* (((v01111./v01TRA111) .* (v00TRA111./v00111) .* (v01001./v01TRA001) .* (v00TRA001./v00001)).^(1/36)) .* (((v01101./v01TRA101) .* (v00TRA101./v00101) .* (v01011./v01TRA011) .* (v00TRA011./v00011)).^(1/72));
# store results in structures
            dv.TR.A.(["r",num2str(r),"_s",num2str(s)]) = dv_TR_A;
            rv.TR.A.(["r",num2str(r),"_s",num2str(s)]) = rv_TR_A;
# calculate the changes in EMPLOYMENT as a weighted average of all TR-related terms - ADDITIVE
            de_TR_A = (e1-e11TRA111)/18 + (e10TRA111-e10111)/18 + (e11001-e11TRA001)/18 + (e10TRA001-e10001)/18 + (e11101-e11TRA101)/36 + (e10TRA101-e10101)/36 + (e11011-e11TRA011)/36 + (e10TRA011-e10011)/36 + (e01110-e01TRA110)/18 + (e00TRA110-e00110)/18 + (e01000-e01TRA000)/18 + (e00TRA000-e0)/18 + (e01100-e01TRA100)/36 + (e00TRA100-e00100)/36 + (e01010-e01TRA010)/36 + (e00TRA010-e00010)/36 + (e11110-e11TRA110)/36 + (e10TRA110-e10110)/36 + (e11000-e11TRA000)/36 + (e10TRA000-e10000)/36 + (e11100-e11TRA100)/72 + (e10TRA100-e10100)/72 + (e11010-e11TRA010)/72 + (e10TRA010-e10010)/72 + (e01111-e01TRA111)/36 + (e00TRA111-e00111)/36 + (e01001-e01TRA001)/36 + (e00TRA001-e00001)/36 + (e01101-e01TRA101)/72 + (e00TRA101-e00101)/72 + (e01011-e01TRA011)/72 + (e00TRA011-e00011)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TR-related terms - MULTIPLICATIVE
            re_TR_A = (((e1./e11TRA111) .* (e10TRA111./e10111) .* (e11001./e11TRA001) .* (e10TRA001./e10001)).^(1/18)) .* (((e11101./e11TRA101) .* (e10TRA101./e10101) .* (e11011./e11TRA011) .* (e10TRA011./e10011)).^(1/36)) .* (((e01110./e01TRA110) .* (e00TRA110./e00110) .* (e01000./e01TRA000) .* (e00TRA000./e0)).^(1/18)) .* (((e01100./e01TRA100) .* (e00TRA100./e00100) .* (e01010./e01TRA010) .* (e00TRA010./e00010)).^(1/36)) .* (((e11110./e11TRA110) .* (e10TRA110./e10110) .* (e11000./e11TRA000) .* (e10TRA000./e10000)).^(1/36)) .* (((e11100./e11TRA100) .* (e10TRA100./e10100) .* (e11010./e11TRA010) .* (e10TRA010./e10010)).^(1/72)) .* (((e01111./e01TRA111) .* (e00TRA111./e00111) .* (e01001./e01TRA001) .* (e00TRA001./e00001)).^(1/36)) .* (((e01101./e01TRA101) .* (e00TRA101./e00101) .* (e01011./e01TRA011) .* (e00TRA011./e00011)).^(1/72));
# store results in structures
            de.TR.A.(["r",num2str(r),"_s",num2str(s)]) = de_TR_A;
            re.TR.A.(["r",num2str(r),"_s",num2str(s)]) = re_TR_A;
# calculate the changes in GDP as a weighted average of all TR-related terms - ADDITIVE
            dy_TR_A = (y11111-y11TRA111)/18 + (y10TRA111-y10111)/18 + (y11001-y11TRA001)/18 + (y10TRA001-y10001)/18 + (y11101-y11TRA101)/36 + (y10TRA101-y10101)/36 + (y11011-y11TRA011)/36 + (y10TRA011-y10011)/36 + (y01110-y01TRA110)/18 + (y00TRA110-y00110)/18 + (y01000-y01TRA000)/18 + (y00TRA000-y00000)/18 + (y01100-y01TRA100)/36 + (y00TRA100-y00100)/36 + (y01010-y01TRA010)/36 + (y00TRA010-y00010)/36 + (y11110-y11TRA110)/36 + (y10TRA110-y10110)/36 + (y11000-y11TRA000)/36 + (y10TRA000-y10000)/36 + (y11100-y11TRA100)/72 + (y10TRA100-y10100)/72 + (y11010-y11TRA010)/72 + (y10TRA010-y10010)/72 + (y01111-y01TRA111)/36 + (y00TRA111-y00111)/36 + (y01001-y01TRA001)/36 + (y00TRA001-y00001)/36 + (y01101-y01TRA101)/72 + (y00TRA101-y00101)/72 + (y01011-y01TRA011)/72 + (y00TRA011-y00011)/72;
# calculate the changes in GDP as a weighted average of all TR-related terms - MULTIPLICATIVE
            ry_TR_A = (((y11111./y11TRA111) .* (y10TRA111./y10111) .* (y11001./y11TRA001) .* (y10TRA001./y10001)).^(1/18)) .* (((y11101./y11TRA101) .* (y10TRA101./y10101) .* (y11011./y11TRA011) .* (y10TRA011./y10011)).^(1/36)) .* (((y01110./y01TRA110) .* (y00TRA110./y00110) .* (y01000./y01TRA000) .* (y00TRA000./y00000)).^(1/18)) .* (((y01100./y01TRA100) .* (y00TRA100./y00100) .* (y01010./y01TRA010) .* (y00TRA010./y00010)).^(1/36)) .* (((y11110./y11TRA110) .* (y10TRA110./y10110) .* (y11000./y11TRA000) .* (y10TRA000./y10000)).^(1/36)) .* (((y11100./y11TRA100) .* (y10TRA100./y10100) .* (y11010./y11TRA010) .* (y10TRA010./y10010)).^(1/72)) .* (((y01111./y01TRA111) .* (y00TRA111./y00111) .* (y01001./y01TRA001) .* (y00TRA001./y00001)).^(1/36)) .* (((y01101./y01TRA101) .* (y00TRA101./y00101) .* (y01011./y01TRA011) .* (y00TRA011./y00011)).^(1/72));
# store results in structures
            dy.TR.A.(["r",num2str(r),"_s",num2str(s)]) = dy_TR_A;
            ry.TR.A.(["r",num2str(r),"_s",num2str(s)]) = ry_TR_A;

# end the if statement for (r != s)
        endif
# end loop for r
    endfor
# clear matrices that are not required until the next loop of s
    clear TCD_A*
# end loop for s
endfor
# clear unnecessary variables to save space (these will be defined again for final demand below, but keep y00000 and y11111)
clear v00* v11* v01* v10* e00* e11* e01* e10* y001* y0001* y00001 y110* y1110* y11110 y01* y10*

# calculate individual terms in the decompositions that do not depend on country pair - FINAL demand
# F111 = F1;
F011 = Fcou0.*(Sk'*Find1*diag(f1));
F100 = Fcou1.*(Sk'*Find0*diag(f0));
# F000 = F0;
F110 = Fcou1.*(Sk'*Find1*diag(f0));
F010 = Fcou0.*(Sk'*Find1*diag(f0));
F101 = Fcou1.*(Sk'*Find0*diag(f1));
F001 = Fcou0.*(Sk'*Find0*diag(f1));
# to simplify calculations, pre-multiply Ti, A and L
Ti1A1L1 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A1),1)))*L1;
Ti0A0L0 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A0),1)))*L0;
Ti1A0L0 = (eye(2340)-diag(sum((1+Ti1).*(Sn'*A0),1)))*L0;
Ti0A1L1 = (eye(2340)-diag(sum((1+Ti0).*(Sn'*A1),1)))*L1;
# to simplify calculations, pre-multiply ec, A and L
ec1L1 = diag(ec1)*L1;
ec0L0 = diag(ec0)*L0;
ec1L0 = diag(ec1)*L0;
ec0L1 = diag(ec0)*L1;

# 32 VALUE ADDED terms for changes in the country of origin of FINAL products
# v***** = v(Tz*,A*,Fcou*,Find,f*)
# v11111 = v1';
v11011 = Ti1A1L1*sum(F011,2);
v11100 = Ti1A1L1*sum(F100,2);
v11000 = Ti1A1L1*sum(F0,2);
v11110 = Ti1A1L1*sum(F110,2);
v11010 = Ti1A1L1*sum(F010,2);
v11101 = Ti1A1L1*sum(F101,2);
v11001 = Ti1A1L1*sum(F001,2);
v00111 = Ti0A0L0*sum(F1,2);
v00011 = Ti0A0L0*sum(F011,2);
v00100 = Ti0A0L0*sum(F100,2);
# v00000 = v0';
v00110 = Ti0A0L0*sum(F110,2);
v00010 = Ti0A0L0*sum(F010,2);
v00101 = Ti0A0L0*sum(F101,2);
v00001 = Ti0A0L0*sum(F001,2);
v10111 = Ti1A0L0*sum(F1,2);
v10011 = Ti1A0L0*sum(F011,2);
v10100 = Ti1A0L0*sum(F100,2);
v10000 = Ti1A0L0*sum(F0,2);
v10110 = Ti1A0L0*sum(F110,2);
v10010 = Ti1A0L0*sum(F010,2);
v10101 = Ti1A0L0*sum(F101,2);
v10001 = Ti1A0L0*sum(F001,2);
v01111 = Ti0A1L1*sum(F1,2);
v01011 = Ti0A1L1*sum(F011,2);
v01100 = Ti0A1L1*sum(F100,2);
v01000 = Ti0A1L1*sum(F0,2);
v01110 = Ti0A1L1*sum(F110,2);
v01010 = Ti0A1L1*sum(F010,2);
v01101 = Ti0A1L1*sum(F101,2);
v01001 = Ti0A1L1*sum(F001,2);

# 32 EMPLOYMENT terms for changes in the country of origin of FINAL products
# e***** = e(ec*,A*,Fcou*,Find,f*)
# e11111 = e1;
e11011 = ec1L1*sum(F011,2);
e11100 = ec1L1*sum(F100,2);
e11000 = ec1L1*sum(F0,2);
e11110 = ec1L1*sum(F110,2);
e11010 = ec1L1*sum(F010,2);
e11101 = ec1L1*sum(F101,2);
e11001 = ec1L1*sum(F001,2);
e00111 = ec0L0*sum(F1,2);
e00011 = ec0L0*sum(F011,2);
e00100 = ec0L0*sum(F100,2);
# e00000 = e0;
e00110 = ec0L0*sum(F110,2);
e00010 = ec0L0*sum(F010,2);
e00101 = ec0L0*sum(F101,2);
e00001 = ec0L0*sum(F001,2);
e10111 = ec1L0*sum(F1,2);
e10011 = ec1L0*sum(F011,2);
e10100 = ec1L0*sum(F100,2);
e10000 = ec1L0*sum(F0,2);
e10110 = ec1L0*sum(F110,2);
e10010 = ec1L0*sum(F010,2);
e10101 = ec1L0*sum(F101,2);
e10001 = ec1L0*sum(F001,2);
e01111 = ec0L1*sum(F1,2);
e01011 = ec0L1*sum(F011,2);
e01100 = ec0L1*sum(F100,2);
e01000 = ec0L1*sum(F0,2);
e01110 = ec0L1*sum(F110,2);
e01010 = ec0L1*sum(F010,2);
e01101 = ec0L1*sum(F101,2);
e01001 = ec0L1*sum(F001,2);

# 32 GDP terms for changes in the country of origin of FINAL products
# calculate (the rest of) taxes less subsidies on final products mf**** = mf(Tf*,Fcou*,Find,f*) for Tf1, Tf0, Fcou*, Find, f*
mf1011 = sum(Tf1.*(Sn'*F011),1)';
mf1100 = sum(Tf1.*(Sn'*F100),1)';
mf1110 = sum(Tf1.*(Sn'*F110),1)';
mf1010 = sum(Tf1.*(Sn'*F010),1)';
mf1101 = sum(Tf1.*(Sn'*F101),1)';
mf1001 = sum(Tf1.*(Sn'*F001),1)';
mf0011 = sum(Tf0.*(Sn'*F011),1)';
mf0100 = sum(Tf0.*(Sn'*F100),1)';
mf0110 = sum(Tf0.*(Sn'*F110),1)';
mf0010 = sum(Tf0.*(Sn'*F010),1)';
mf0101 = sum(Tf0.*(Sn'*F101),1)';
mf0001 = sum(Tf0.*(Sn'*F001),1)';
# to simplify calculations, pre-multiply A and L
A1L1 = (eye(2340)-diag(sum(A1,1)))*L1;
A0L0 = (eye(2340)-diag(sum(A0,1)))*L0;
# y***** = y(Tf*,A*,Fcou*,Find,f*)
# y11111 = Sn'*(x1 - sum(Z1,1))' + mf1111;
y11011 = Sn'*A1L1*sum(F011,2) + mf1011;
y11100 = Sn'*A1L1*sum(F100,2) + mf1100;
y11000 = Sn'*A1L1*sum(F0,2) + mf1000;
y11110 = Sn'*A1L1*sum(F110,2) + mf1110;
y11010 = Sn'*A1L1*sum(F010,2) + mf1010;
y11101 = Sn'*A1L1*sum(F101,2) + mf1101;
y11001 = Sn'*A1L1*sum(F001,2) + mf1001;
y00111 = Sn'*A0L0*sum(F1,2) + mf0111;
y00011 = Sn'*A0L0*sum(F011,2) + mf0011;
y00100 = Sn'*A0L0*sum(F100,2) + mf0100;
# y00000 = Sn'*(x0 - sum(Z0,1))' + mf0000;
y00110 = Sn'*A0L0*sum(F110,2) + mf0110;
y00010 = Sn'*A0L0*sum(F010,2) + mf0010;
y00101 = Sn'*A0L0*sum(F101,2) + mf0101;
y00001 = Sn'*A0L0*sum(F001,2) + mf0001;
y10111 = Sn'*A0L0*sum(F1,2) + mf1111;
y10011 = Sn'*A0L0*sum(F011,2) + mf1011;
y10100 = Sn'*A0L0*sum(F100,2) + mf1100;
y10000 = Sn'*A0L0*sum(F0,2) + mf1000;
y10110 = Sn'*A0L0*sum(F110,2) + mf1110;
y10010 = Sn'*A0L0*sum(F010,2) + mf1010;
y10101 = Sn'*A0L0*sum(F101,2) + mf1101;
y10001 = Sn'*A0L0*sum(F001,2) + mf1001;
y01111 = Sn'*A1L1*sum(F1,2) + mf0111;
y01011 = Sn'*A1L1*sum(F011,2) + mf0011;
y01100 = Sn'*A1L1*sum(F100,2) + mf0100;
y01000 = Sn'*A1L1*sum(F0,2) + mf0000;
y01110 = Sn'*A1L1*sum(F110,2) + mf0110;
y01010 = Sn'*A1L1*sum(F010,2) + mf0010;
y01101 = Sn'*A1L1*sum(F101,2) + mf0101;
y01001 = Sn'*A1L1*sum(F001,2) + mf0001;

# clear unnecessary F*** terms to save space
clear F01* F10*

# start again the loop for s, in the OECD ICIO 2018: Korea is No.19 and USA is No.36
for s = [19,36]

# run a bi-proportional balancing (equivalent to two RAS runs) to a matrix of ones – final products
# r - partner (country of origin), s - importer (country in focus), i - industry of origin, j - producer (industry in focus)
# for each sth column in dFcou_neg and dFcou_pos, generate a block-diagonal matrix with N K x K matrices in its diagonal blocks; this implicitly involves transition from the country-industry to the industry-country classification, so each block in the matrix corresponds to product i
    TCD_F0 = kron(speye(36),ones(65));
# to apply the bi-proportional method of balancing, I adjust the block-diagonal matrix of ones TCD_F0 by scaling it row-wise and column-wise as in he first and second RAS runs (only two RAS runs are enough in this case!), so that it conforms with the row and column totals
# calculate a column vector of row multipliers
    rRAS = (P1*dFcou_pos(:,s))./sum(TCD_F0, 2);
# scale row-wise
    TCD_F1 = diag(rRAS)*TCD_F0;
# calculate a row vector of column multipliers
    sRAS = ((-P1*dFcou_neg(:,s))')./sum(TCD_F1, 1);
# address possible "NaN" in sRAS, replacing those with 0
	sRAS(isnan(sRAS)==1) = 0;
# scale column-wise
    TCD_F1 = TCD_F1*diag(sRAS);
# TCD_F1 is a block-diagonal matrix that shows the re-distribution between countries of origin (r) of each final product (i) used in the selected country (s)
# reorganize TCD_F1 into one block column and stack to finally obtain a KN x K matrix
    TCD_F = sparse(2340, 65);
    for n = 1:36
        TCD_F(1+65*(n-1):65+65*(n-1),:) = TCD_F1(1+65*(n-1):65+65*(n-1),1+65*(n-1):65+65*(n-1));
    endfor
    clear TCD_F1

# extract the information from TCD_F on trade creation (diversion, contraction) between country r and country s
# reshuffle rows and columns with the permutation matrix P, revert to the original country-industry classification
    TCD_F = P1'*TCD_F;
    for r = [36,19]
# ignore a case where r = s
        if (r != s)
# construct a matrix of trade creation between countries s and r
            TC_F = sparse(2340, 65);
# copy the r-s block from TCD_F to the r-s block of TC_F
            TC_F(1+36*(r-1):36+36*(r-1),s) = TCD_F(1+36*(r-1):36+36*(r-1),s);
# copy the r-s block of TCD_A to the s-s block of TC_A, changing sign
            TC_F(1+36*(s-1):36+36*(s-1),s) = -TCD_F(1+36*(r-1):36+36*(r-1),s);
# construct a matrix of trade diversion in between third countries and r
            TDin_F = sparse(2340, 65);
# block-transpose the row of r blocks from TCD_F, changing sign
            for t = 1:65
                TDin_F(1+36*(t-1):36+36*(t-1),s) = -TCD_F(1+36*(r-1):36+36*(r-1),t);
            endfor
# calculate the TDin for country r as aggregate TDin in all r-t pairs less TC in the r-s pair
            TDin_F(1+36*(r-1):36+36*(r-1),s) = sum(TCD_F(1+36*(r-1):36+36*(r-1),:),2) - TCD_F(1+36*(r-1):36+36*(r-1),s);
# set the s-s block to zero
            TDin_F(1+36*(s-1):36+36*(s-1),s) = zeros(36, 1);
# construct a matrix of trade diversion out between third countries and r
            TDout_F = sparse(2340, 65);
# copy the column of r from TCD_F to the column of s in TDout_F
            TDout_F(:,s) = TCD_F(:,r);
# calculate TDout for country r as aggregate TDout in all t-r pairs less TR in the r-s pair, and change sign
            TDout_F(1+36*(r-1):36+36*(r-1),s) = -(Sk*TCD_F(:,r) - TCD_F(1+36*(s-1):36+36*(s-1),r));
# set the s-s block to zero
            TDout_F(1+36*(s-1):36+36*(s-1),s) = zeros(36, 1);
# construct a matrix of trade contraction (import substitution) between s and r
            TR_F = sparse(2340, 65);
# copy the s-r block from TCD_F to the s-s block of TR_F
            TR_F(1+36*(s-1):36+36*(s-1),s) = TCD_F(1+36*(s-1):36+36*(s-1),r);
# copy the s-r block from TCD_F to the r-s block of TR_F, changing sign
            TR_F(1+36*(r-1):36+36*(r-1),s) = -TCD_F(1+36*(s-1):36+36*(s-1),r);

# the rest of the (bilateral) component mattices for SDA - TC
            F1TCF11 = (Fcou1-TC_F).*(Sk'*Find1*diag(f1));
            F0TCF11 = (Fcou0+TC_F).*(Sk'*Find1*diag(f1));
            F1TCF00 = (Fcou1-TC_F).*(Sk'*Find0*diag(f0));
            F0TCF00 = (Fcou0+TC_F).*(Sk'*Find0*diag(f0));
            F1TCF10 = (Fcou1-TC_F).*(Sk'*Find1*diag(f0));
            F0TCF10 = (Fcou0+TC_F).*(Sk'*Find1*diag(f0));
            F1TCF01 = (Fcou1-TC_F).*(Sk'*Find0*diag(f1));
            F0TCF01 = (Fcou0+TC_F).*(Sk'*Find0*diag(f1));
# the rest 32 VALUE ADDED terms for changes in the country of orgin of FINAL demand
# v***TCF** = v(Tz*,A*,Fcou*,Find*,f*)
            v111TCF11 = Ti1A1L1*sum(F1TCF11,2);
            v110TCF11 = Ti1A1L1*sum(F0TCF11,2);
            v111TCF00 = Ti1A1L1*sum(F1TCF00,2);
            v110TCF00 = Ti1A1L1*sum(F0TCF00,2);
            v111TCF10 = Ti1A1L1*sum(F1TCF10,2);
            v110TCF10 = Ti1A1L1*sum(F0TCF10,2);
            v111TCF01 = Ti1A1L1*sum(F1TCF01,2);
            v110TCF01 = Ti1A1L1*sum(F0TCF01,2);
            v001TCF11 = Ti0A0L0*sum(F1TCF11,2);
            v000TCF11 = Ti0A0L0*sum(F0TCF11,2);
            v001TCF00 = Ti0A0L0*sum(F1TCF00,2);
            v000TCF00 = Ti0A0L0*sum(F0TCF00,2);
            v001TCF10 = Ti0A0L0*sum(F1TCF10,2);
            v000TCF10 = Ti0A0L0*sum(F0TCF10,2);
            v001TCF01 = Ti0A0L0*sum(F1TCF01,2);
            v000TCF01 = Ti0A0L0*sum(F0TCF01,2);
            v101TCF11 = Ti1A0L0*sum(F1TCF11,2);
            v100TCF11 = Ti1A0L0*sum(F0TCF11,2);
            v101TCF00 = Ti1A0L0*sum(F1TCF00,2);
            v100TCF00 = Ti1A0L0*sum(F0TCF00,2);
            v101TCF10 = Ti1A0L0*sum(F1TCF10,2);
            v100TCF10 = Ti1A0L0*sum(F0TCF10,2);
            v101TCF01 = Ti1A0L0*sum(F1TCF01,2);
            v100TCF01 = Ti1A0L0*sum(F0TCF01,2);
            v011TCF11 = Ti0A1L1*sum(F1TCF11,2);
            v010TCF11 = Ti0A1L1*sum(F0TCF11,2);
            v011TCF00 = Ti0A1L1*sum(F1TCF00,2);
            v010TCF00 = Ti0A1L1*sum(F0TCF00,2);
            v011TCF10 = Ti0A1L1*sum(F1TCF10,2);
            v010TCF10 = Ti0A1L1*sum(F0TCF10,2);
            v011TCF01 = Ti0A1L1*sum(F1TCF01,2);
            v010TCF01 = Ti0A1L1*sum(F0TCF01,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of FINAL demand
# e***TCF** = e(ec*,A*,Fcou*,Find*,f*)
            e111TCF11 = ec1L1*sum(F1TCF11,2);
            e110TCF11 = ec1L1*sum(F0TCF11,2);
            e111TCF00 = ec1L1*sum(F1TCF00,2);
            e110TCF00 = ec1L1*sum(F0TCF00,2);
            e111TCF10 = ec1L1*sum(F1TCF10,2);
            e110TCF10 = ec1L1*sum(F0TCF10,2);
            e111TCF01 = ec1L1*sum(F1TCF01,2);
            e110TCF01 = ec1L1*sum(F0TCF01,2);
            e001TCF11 = ec0L0*sum(F1TCF11,2);
            e000TCF11 = ec0L0*sum(F0TCF11,2);
            e001TCF00 = ec0L0*sum(F1TCF00,2);
            e000TCF00 = ec0L0*sum(F0TCF00,2);
            e001TCF10 = ec0L0*sum(F1TCF10,2);
            e000TCF10 = ec0L0*sum(F0TCF10,2);
            e001TCF01 = ec0L0*sum(F1TCF01,2);
            e000TCF01 = ec0L0*sum(F0TCF01,2);
            e101TCF11 = ec1L0*sum(F1TCF11,2);
            e100TCF11 = ec1L0*sum(F0TCF11,2);
            e101TCF00 = ec1L0*sum(F1TCF00,2);
            e100TCF00 = ec1L0*sum(F0TCF00,2);
            e101TCF10 = ec1L0*sum(F1TCF10,2);
            e100TCF10 = ec1L0*sum(F0TCF10,2);
            e101TCF01 = ec1L0*sum(F1TCF01,2);
            e100TCF01 = ec1L0*sum(F0TCF01,2);
            e011TCF11 = ec0L1*sum(F1TCF11,2);
            e010TCF11 = ec0L1*sum(F0TCF11,2);
            e011TCF00 = ec0L1*sum(F1TCF00,2);
            e010TCF00 = ec0L1*sum(F0TCF00,2);
            e011TCF10 = ec0L1*sum(F1TCF10,2);
            e010TCF10 = ec0L1*sum(F0TCF10,2);
            e011TCF01 = ec0L1*sum(F1TCF01,2);
            e010TCF01 = ec0L1*sum(F0TCF01,2);
# the rest 16 TAXES LESS SUBSIDIES terms for changes in the country of origin of FINAL products
# mf**TCF** = mf(Tf*,Fcou*,Find*,f*)
            mf11TCF11 = sum(Tf1.*(Sn'*F1TCF11),1)';
            mf10TCF11 = sum(Tf1.*(Sn'*F0TCF11),1)';
            mf11TCF00 = sum(Tf1.*(Sn'*F1TCF00),1)';
            mf10TCF00 = sum(Tf1.*(Sn'*F0TCF00),1)';
            mf11TCF10 = sum(Tf1.*(Sn'*F1TCF10),1)';
            mf10TCF10 = sum(Tf1.*(Sn'*F0TCF10),1)';
            mf11TCF01 = sum(Tf1.*(Sn'*F1TCF01),1)';
            mf10TCF01 = sum(Tf1.*(Sn'*F0TCF01),1)';
            mf01TCF11 = sum(Tf0.*(Sn'*F1TCF11),1)';
            mf00TCF11 = sum(Tf0.*(Sn'*F0TCF11),1)';
            mf01TCF00 = sum(Tf0.*(Sn'*F1TCF00),1)';
            mf00TCF00 = sum(Tf0.*(Sn'*F0TCF00),1)';
            mf01TCF10 = sum(Tf0.*(Sn'*F1TCF10),1)';
            mf00TCF10 = sum(Tf0.*(Sn'*F0TCF10),1)';
            mf01TCF01 = sum(Tf0.*(Sn'*F1TCF01),1)';
            mf00TCF01 = sum(Tf0.*(Sn'*F0TCF01),1)';
# the rest 32 GDP terms for changes in the country of origin of FINAL products
# y***TCF** = y(Tf*,A,Fcou*,Find*,f*)
            y111TCF11 = Sn'*A1L1*sum(F1TCF11,2) + mf11TCF11;
            y110TCF11 = Sn'*A1L1*sum(F0TCF11,2) + mf10TCF11;
            y111TCF00 = Sn'*A1L1*sum(F1TCF00,2) + mf11TCF00;
            y110TCF00 = Sn'*A1L1*sum(F0TCF00,2) + mf10TCF00;
            y111TCF10 = Sn'*A1L1*sum(F1TCF10,2) + mf11TCF10;
            y110TCF10 = Sn'*A1L1*sum(F0TCF10,2) + mf10TCF10;
            y111TCF01 = Sn'*A1L1*sum(F1TCF01,2) + mf11TCF01;
            y110TCF01 = Sn'*A1L1*sum(F0TCF01,2) + mf10TCF01;
            y001TCF11 = Sn'*A0L0*sum(F1TCF11,2) + mf01TCF11;
            y000TCF11 = Sn'*A0L0*sum(F0TCF11,2) + mf00TCF11;
            y001TCF00 = Sn'*A0L0*sum(F1TCF00,2) + mf01TCF00;
            y000TCF00 = Sn'*A0L0*sum(F0TCF00,2) + mf00TCF00;
            y001TCF10 = Sn'*A0L0*sum(F1TCF10,2) + mf01TCF10;
            y000TCF10 = Sn'*A0L0*sum(F0TCF10,2) + mf00TCF10;
            y001TCF01 = Sn'*A0L0*sum(F1TCF01,2) + mf01TCF01;
            y000TCF01 = Sn'*A0L0*sum(F0TCF01,2) + mf00TCF01;
            y101TCF11 = Sn'*A0L0*sum(F1TCF11,2) + mf11TCF11;
            y100TCF11 = Sn'*A0L0*sum(F0TCF11,2) + mf10TCF11;
            y101TCF00 = Sn'*A0L0*sum(F1TCF00,2) + mf11TCF00;
            y100TCF00 = Sn'*A0L0*sum(F0TCF00,2) + mf10TCF00;
            y101TCF10 = Sn'*A0L0*sum(F1TCF10,2) + mf11TCF10;
            y100TCF10 = Sn'*A0L0*sum(F0TCF10,2) + mf10TCF10;
            y101TCF01 = Sn'*A0L0*sum(F1TCF01,2) + mf11TCF01;
            y100TCF01 = Sn'*A0L0*sum(F0TCF01,2) + mf10TCF01;
            y011TCF11 = Sn'*A1L1*sum(F1TCF11,2) + mf01TCF11;
            y010TCF11 = Sn'*A1L1*sum(F0TCF11,2) + mf00TCF11;
            y011TCF00 = Sn'*A1L1*sum(F1TCF00,2) + mf01TCF00;
            y010TCF00 = Sn'*A1L1*sum(F0TCF00,2) + mf00TCF00;
            y011TCF10 = Sn'*A1L1*sum(F1TCF10,2) + mf01TCF10;
            y010TCF10 = Sn'*A1L1*sum(F0TCF10,2) + mf00TCF10;
            y011TCF01 = Sn'*A1L1*sum(F1TCF01,2) + mf01TCF01;
            y010TCF01 = Sn'*A1L1*sum(F0TCF01,2) + mf00TCF01;
# clear matrices that are not required until the next loop of r
            clear TC_F F0TCF* F1TCF*
# calculate the changes in VALUE ADDED as a weighted average of all TC-related terms - ADDITIVE
            dv_TC_F = (v1'-v111TCF11)/18 + (v110TCF11-v11011)/18 + (v11100-v111TCF00)/18 + (v110TCF00-v11000)/18 + (v11110-v111TCF10)/36 + (v110TCF10-v11010)/36 + (v11101-v111TCF01)/36 + (v110TCF01-v11001)/36 + (v00111-v001TCF11)/18 + (v000TCF11-v00011)/18 + (v00100-v001TCF00)/18 + (v000TCF00-v0')/18 + (v00110-v001TCF10)/36 + (v000TCF10-v00010)/36 + (v00101-v001TCF01)/36 + (v000TCF01-v00001)/36 + (v10111-v101TCF11)/36 + (v100TCF11-v10011)/36 + (v10100-v101TCF00)/36 + (v100TCF00-v10000)/36 + (v10110-v101TCF10)/72 + (v100TCF10-v10010)/72 + (v10101-v101TCF01)/72 + (v100TCF01-v10001)/72 + (v01111-v011TCF11)/36 + (v010TCF11-v01011)/36 + (v01100-v011TCF00)/36 + (v010TCF00-v01000)/36 + (v01110-v011TCF10)/72 + (v010TCF10-v01010)/72 + (v01101-v011TCF01)/72 + (v010TCF01-v01001)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TC-related terms - MULTIPLICATIVE
            rv_TC_F = (((v1'./v111TCF11) .* (v110TCF11./v11011) .* (v11100./v111TCF00) .* (v110TCF00./v11000)).^(1/18)) .* (((v11110./v111TCF10) .* (v110TCF10./v11010) .* (v11101./v111TCF01) .* (v110TCF01./v11001)).^(1/36)) .* (((v00111./v001TCF11) .* (v000TCF11./v00011) .* (v00100./v001TCF00) .* (v000TCF00./v0')).^(1/18)) .* (((v00110./v001TCF10) .* (v000TCF10./v00010) .* (v00101./v001TCF01) .* (v000TCF01./v00001)).^(1/36)) .* (((v10111./v101TCF11) .* (v100TCF11./v10011) .* (v10100./v101TCF00) .* (v100TCF00./v10000)).^(1/36)) .* (((v10110./v101TCF10) .* (v100TCF10./v10010) .* (v10101./v101TCF01) .* (v100TCF01./v10001)).^(1/72)) .* (((v01111./v011TCF11) .* (v010TCF11./v01011) .* (v01100./v011TCF00) .* (v010TCF00./v01000)).^(1/36)) .* (((v01110./v011TCF10) .* (v010TCF10./v01010) .* (v01101./v011TCF01) .* (v010TCF01./v01001)).^(1/72));
# store results in structures
            dv.TC.F.(["r",num2str(r),"_s",num2str(s)]) = dv_TC_F;
            rv.TC.F.(["r",num2str(r),"_s",num2str(s)]) = rv_TC_F;
# calculate the changes in EMPLOYMENT as a weighted average of all TC-related terms - ADDITIVE
            de_TC_F = (e1-e111TCF11)/18 + (e110TCF11-e11011)/18 + (e11100-e111TCF00)/18 + (e110TCF00-e11000)/18 + (e11110-e111TCF10)/36 + (e110TCF10-e11010)/36 + (e11101-e111TCF01)/36 + (e110TCF01-e11001)/36 + (e00111-e001TCF11)/18 + (e000TCF11-e00011)/18 + (e00100-e001TCF00)/18 + (e000TCF00-e0)/18 + (e00110-e001TCF10)/36 + (e000TCF10-e00010)/36 + (e00101-e001TCF01)/36 + (e000TCF01-e00001)/36 + (e10111-e101TCF11)/36 + (e100TCF11-e10011)/36 + (e10100-e101TCF00)/36 + (e100TCF00-e10000)/36 + (e10110-e101TCF10)/72 + (e100TCF10-e10010)/72 + (e10101-e101TCF01)/72 + (e100TCF01-e10001)/72 + (e01111-e011TCF11)/36 + (e010TCF11-e01011)/36 + (e01100-e011TCF00)/36 + (e010TCF00-e01000)/36 + (e01110-e011TCF10)/72 + (e010TCF10-e01010)/72 + (e01101-e011TCF01)/72 + (e010TCF01-e01001)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TC-related terms - MULTIPLICATIVE
            re_TC_F = (((e1./e111TCF11) .* (e110TCF11./e11011) .* (e11100./e111TCF00) .* (e110TCF00./e11000)).^(1/18)) .* (((e11110./e111TCF10) .* (e110TCF10./e11010) .* (e11101./e111TCF01) .* (e110TCF01./e11001)).^(1/36)) .* (((e00111./e001TCF11) .* (e000TCF11./e00011) .* (e00100./e001TCF00) .* (e000TCF00./e0)).^(1/18)) .* (((e00110./e001TCF10) .* (e000TCF10./e00010) .* (e00101./e001TCF01) .* (e000TCF01./e00001)).^(1/36)) .* (((e10111./e101TCF11) .* (e100TCF11./e10011) .* (e10100./e101TCF00) .* (e100TCF00./e10000)).^(1/36)) .* (((e10110./e101TCF10) .* (e100TCF10./e10010) .* (e10101./e101TCF01) .* (e100TCF01./e10001)).^(1/72)) .* (((e01111./e011TCF11) .* (e010TCF11./e01011) .* (e01100./e011TCF00) .* (e010TCF00./e01000)).^(1/36)) .* (((e01110./e011TCF10) .* (e010TCF10./e01010) .* (e01101./e011TCF01) .* (e010TCF01./e01001)).^(1/72));
# store results in structures
            de.TC.F.(["r",num2str(r),"_s",num2str(s)]) = de_TC_F;
            re.TC.F.(["r",num2str(r),"_s",num2str(s)]) = re_TC_F;
# calculate the changes in GDP as a weighted average of all TC-related terms - ADDITIVE
            dy_TC_F = (y11111-y111TCF11)/18 + (y110TCF11-y11011)/18 + (y11100-y111TCF00)/18 + (y110TCF00-y11000)/18 + (y11110-y111TCF10)/36 + (y110TCF10-y11010)/36 + (y11101-y111TCF01)/36 + (y110TCF01-y11001)/36 + (y00111-y001TCF11)/18 + (y000TCF11-y00011)/18 + (y00100-y001TCF00)/18 + (y000TCF00-y00000)/18 + (y00110-y001TCF10)/36 + (y000TCF10-y00010)/36 + (y00101-y001TCF01)/36 + (y000TCF01-y00001)/36 + (y10111-y101TCF11)/36 + (y100TCF11-y10011)/36 + (y10100-y101TCF00)/36 + (y100TCF00-y10000)/36 + (y10110-y101TCF10)/72 + (y100TCF10-y10010)/72 + (y10101-y101TCF01)/72 + (y100TCF01-y10001)/72 + (y01111-y011TCF11)/36 + (y010TCF11-y01011)/36 + (y01100-y011TCF00)/36 + (y010TCF00-y01000)/36 + (y01110-y011TCF10)/72 + (y010TCF10-y01010)/72 + (y01101-y011TCF01)/72 + (y010TCF01-y01001)/72;
# calculate the changes in GDP as a weighted average of all TC-related terms - MULTIPLICATIVE
            ry_TC_F = (((y11111./y111TCF11) .* (y110TCF11./y11011) .* (y11100./y111TCF00) .* (y110TCF00./y11000)).^(1/18)) .* (((y11110./y111TCF10) .* (y110TCF10./y11010) .* (y11101./y111TCF01) .* (y110TCF01./y11001)).^(1/36)) .* (((y00111./y001TCF11) .* (y000TCF11./y00011) .* (y00100./y001TCF00) .* (y000TCF00./y00000)).^(1/18)) .* (((y00110./y001TCF10) .* (y000TCF10./y00010) .* (y00101./y001TCF01) .* (y000TCF01./y00001)).^(1/36)) .* (((y10111./y101TCF11) .* (y100TCF11./y10011) .* (y10100./y101TCF00) .* (y100TCF00./y10000)).^(1/36)) .* (((y10110./y101TCF10) .* (y100TCF10./y10010) .* (y10101./y101TCF01) .* (y100TCF01./y10001)).^(1/72)) .* (((y01111./y011TCF11) .* (y010TCF11./y01011) .* (y01100./y011TCF00) .* (y010TCF00./y01000)).^(1/36)) .* (((y01110./y011TCF10) .* (y010TCF10./y01010) .* (y01101./y011TCF01) .* (y010TCF01./y01001)).^(1/72));
# store results in structures
            dy.TC.F.(["r",num2str(r),"_s",num2str(s)]) = dy_TC_F;
            ry.TC.F.(["r",num2str(r),"_s",num2str(s)]) = ry_TC_F;

# the rest of the (bilateral) component matrices - TDin
            F1TDinF11 = (Fcou1-TDin_F).*(Sk'*Find1*diag(f1));
            F0TDinF11 = (Fcou0+TDin_F).*(Sk'*Find1*diag(f1));
            F1TDinF00 = (Fcou1-TDin_F).*(Sk'*Find0*diag(f0));
            F0TDinF00 = (Fcou0+TDin_F).*(Sk'*Find0*diag(f0));
            F1TDinF10 = (Fcou1-TDin_F).*(Sk'*Find1*diag(f0));
            F0TDinF10 = (Fcou0+TDin_F).*(Sk'*Find1*diag(f0));
            F1TDinF01 = (Fcou1-TDin_F).*(Sk'*Find0*diag(f1));
            F0TDinF01 = (Fcou0+TDin_F).*(Sk'*Find0*diag(f1));
# the rest 32 VALUE ADDED terms for changes in the country of orgin of FINAL demand
# v***TDinF** = v(Tz*,A*,Fcou*,Find*,f*)
            v111TDinF11 = Ti1A1L1*sum(F1TDinF11,2);
            v110TDinF11 = Ti1A1L1*sum(F0TDinF11,2);
            v111TDinF00 = Ti1A1L1*sum(F1TDinF00,2);
            v110TDinF00 = Ti1A1L1*sum(F0TDinF00,2);
            v111TDinF10 = Ti1A1L1*sum(F1TDinF10,2);
            v110TDinF10 = Ti1A1L1*sum(F0TDinF10,2);
            v111TDinF01 = Ti1A1L1*sum(F1TDinF01,2);
            v110TDinF01 = Ti1A1L1*sum(F0TDinF01,2);
            v001TDinF11 = Ti0A0L0*sum(F1TDinF11,2);
            v000TDinF11 = Ti0A0L0*sum(F0TDinF11,2);
            v001TDinF00 = Ti0A0L0*sum(F1TDinF00,2);
            v000TDinF00 = Ti0A0L0*sum(F0TDinF00,2);
            v001TDinF10 = Ti0A0L0*sum(F1TDinF10,2);
            v000TDinF10 = Ti0A0L0*sum(F0TDinF10,2);
            v001TDinF01 = Ti0A0L0*sum(F1TDinF01,2);
            v000TDinF01 = Ti0A0L0*sum(F0TDinF01,2);
            v101TDinF11 = Ti1A0L0*sum(F1TDinF11,2);
            v100TDinF11 = Ti1A0L0*sum(F0TDinF11,2);
            v101TDinF00 = Ti1A0L0*sum(F1TDinF00,2);
            v100TDinF00 = Ti1A0L0*sum(F0TDinF00,2);
            v101TDinF10 = Ti1A0L0*sum(F1TDinF10,2);
            v100TDinF10 = Ti1A0L0*sum(F0TDinF10,2);
            v101TDinF01 = Ti1A0L0*sum(F1TDinF01,2);
            v100TDinF01 = Ti1A0L0*sum(F0TDinF01,2);
            v011TDinF11 = Ti0A1L1*sum(F1TDinF11,2);
            v010TDinF11 = Ti0A1L1*sum(F0TDinF11,2);
            v011TDinF00 = Ti0A1L1*sum(F1TDinF00,2);
            v010TDinF00 = Ti0A1L1*sum(F0TDinF00,2);
            v011TDinF10 = Ti0A1L1*sum(F1TDinF10,2);
            v010TDinF10 = Ti0A1L1*sum(F0TDinF10,2);
            v011TDinF01 = Ti0A1L1*sum(F1TDinF01,2);
            v010TDinF01 = Ti0A1L1*sum(F0TDinF01,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of FINAL demand
# e***TDinF** = e(ec*,A*,Fcou*,Find*,f*)
            e111TDinF11 = ec1L1*sum(F1TDinF11,2);
            e110TDinF11 = ec1L1*sum(F0TDinF11,2);
            e111TDinF00 = ec1L1*sum(F1TDinF00,2);
            e110TDinF00 = ec1L1*sum(F0TDinF00,2);
            e111TDinF10 = ec1L1*sum(F1TDinF10,2);
            e110TDinF10 = ec1L1*sum(F0TDinF10,2);
            e111TDinF01 = ec1L1*sum(F1TDinF01,2);
            e110TDinF01 = ec1L1*sum(F0TDinF01,2);
            e001TDinF11 = ec0L0*sum(F1TDinF11,2);
            e000TDinF11 = ec0L0*sum(F0TDinF11,2);
            e001TDinF00 = ec0L0*sum(F1TDinF00,2);
            e000TDinF00 = ec0L0*sum(F0TDinF00,2);
            e001TDinF10 = ec0L0*sum(F1TDinF10,2);
            e000TDinF10 = ec0L0*sum(F0TDinF10,2);
            e001TDinF01 = ec0L0*sum(F1TDinF01,2);
            e000TDinF01 = ec0L0*sum(F0TDinF01,2);
            e101TDinF11 = ec1L0*sum(F1TDinF11,2);
            e100TDinF11 = ec1L0*sum(F0TDinF11,2);
            e101TDinF00 = ec1L0*sum(F1TDinF00,2);
            e100TDinF00 = ec1L0*sum(F0TDinF00,2);
            e101TDinF10 = ec1L0*sum(F1TDinF10,2);
            e100TDinF10 = ec1L0*sum(F0TDinF10,2);
            e101TDinF01 = ec1L0*sum(F1TDinF01,2);
            e100TDinF01 = ec1L0*sum(F0TDinF01,2);
            e011TDinF11 = ec0L1*sum(F1TDinF11,2);
            e010TDinF11 = ec0L1*sum(F0TDinF11,2);
            e011TDinF00 = ec0L1*sum(F1TDinF00,2);
            e010TDinF00 = ec0L1*sum(F0TDinF00,2);
            e011TDinF10 = ec0L1*sum(F1TDinF10,2);
            e010TDinF10 = ec0L1*sum(F0TDinF10,2);
            e011TDinF01 = ec0L1*sum(F1TDinF01,2);
            e010TDinF01 = ec0L1*sum(F0TDinF01,2);
# the rest 16 TAXES LESS SUBSIDIES terms for changes in the country of origin of FINAL products
# mf**TDinF** = mf(Tf*,Fcou*,Find*,f*)
            mf11TDinF11 = sum(Tf1.*(Sn'*F1TDinF11),1)';
            mf10TDinF11 = sum(Tf1.*(Sn'*F0TDinF11),1)';
            mf11TDinF00 = sum(Tf1.*(Sn'*F1TDinF00),1)';
            mf10TDinF00 = sum(Tf1.*(Sn'*F0TDinF00),1)';
            mf11TDinF10 = sum(Tf1.*(Sn'*F1TDinF10),1)';
            mf10TDinF10 = sum(Tf1.*(Sn'*F0TDinF10),1)';
            mf11TDinF01 = sum(Tf1.*(Sn'*F1TDinF01),1)';
            mf10TDinF01 = sum(Tf1.*(Sn'*F0TDinF01),1)';
            mf01TDinF11 = sum(Tf0.*(Sn'*F1TDinF11),1)';
            mf00TDinF11 = sum(Tf0.*(Sn'*F0TDinF11),1)';
            mf01TDinF00 = sum(Tf0.*(Sn'*F1TDinF00),1)';
            mf00TDinF00 = sum(Tf0.*(Sn'*F0TDinF00),1)';
            mf01TDinF10 = sum(Tf0.*(Sn'*F1TDinF10),1)';
            mf00TDinF10 = sum(Tf0.*(Sn'*F0TDinF10),1)';
            mf01TDinF01 = sum(Tf0.*(Sn'*F1TDinF01),1)';
            mf00TDinF01 = sum(Tf0.*(Sn'*F0TDinF01),1)';
# the rest 32 GDP terms for changes in the country of origin of FINAL products
# y***TDinF** = y(Tf*,A,Fcou*,Find*,f*)
            y111TDinF11 = Sn'*A1L1*sum(F1TDinF11,2) + mf11TDinF11;
            y110TDinF11 = Sn'*A1L1*sum(F0TDinF11,2) + mf10TDinF11;
            y111TDinF00 = Sn'*A1L1*sum(F1TDinF00,2) + mf11TDinF00;
            y110TDinF00 = Sn'*A1L1*sum(F0TDinF00,2) + mf10TDinF00;
            y111TDinF10 = Sn'*A1L1*sum(F1TDinF10,2) + mf11TDinF10;
            y110TDinF10 = Sn'*A1L1*sum(F0TDinF10,2) + mf10TDinF10;
            y111TDinF01 = Sn'*A1L1*sum(F1TDinF01,2) + mf11TDinF01;
            y110TDinF01 = Sn'*A1L1*sum(F0TDinF01,2) + mf10TDinF01;
            y001TDinF11 = Sn'*A0L0*sum(F1TDinF11,2) + mf01TDinF11;
            y000TDinF11 = Sn'*A0L0*sum(F0TDinF11,2) + mf00TDinF11;
            y001TDinF00 = Sn'*A0L0*sum(F1TDinF00,2) + mf01TDinF00;
            y000TDinF00 = Sn'*A0L0*sum(F0TDinF00,2) + mf00TDinF00;
            y001TDinF10 = Sn'*A0L0*sum(F1TDinF10,2) + mf01TDinF10;
            y000TDinF10 = Sn'*A0L0*sum(F0TDinF10,2) + mf00TDinF10;
            y001TDinF01 = Sn'*A0L0*sum(F1TDinF01,2) + mf01TDinF01;
            y000TDinF01 = Sn'*A0L0*sum(F0TDinF01,2) + mf00TDinF01;
            y101TDinF11 = Sn'*A0L0*sum(F1TDinF11,2) + mf11TDinF11;
            y100TDinF11 = Sn'*A0L0*sum(F0TDinF11,2) + mf10TDinF11;
            y101TDinF00 = Sn'*A0L0*sum(F1TDinF00,2) + mf11TDinF00;
            y100TDinF00 = Sn'*A0L0*sum(F0TDinF00,2) + mf10TDinF00;
            y101TDinF10 = Sn'*A0L0*sum(F1TDinF10,2) + mf11TDinF10;
            y100TDinF10 = Sn'*A0L0*sum(F0TDinF10,2) + mf10TDinF10;
            y101TDinF01 = Sn'*A0L0*sum(F1TDinF01,2) + mf11TDinF01;
            y100TDinF01 = Sn'*A0L0*sum(F0TDinF01,2) + mf10TDinF01;
            y011TDinF11 = Sn'*A1L1*sum(F1TDinF11,2) + mf01TDinF11;
            y010TDinF11 = Sn'*A1L1*sum(F0TDinF11,2) + mf00TDinF11;
            y011TDinF00 = Sn'*A1L1*sum(F1TDinF00,2) + mf01TDinF00;
            y010TDinF00 = Sn'*A1L1*sum(F0TDinF00,2) + mf00TDinF00;
            y011TDinF10 = Sn'*A1L1*sum(F1TDinF10,2) + mf01TDinF10;
            y010TDinF10 = Sn'*A1L1*sum(F0TDinF10,2) + mf00TDinF10;
            y011TDinF01 = Sn'*A1L1*sum(F1TDinF01,2) + mf01TDinF01;
            y010TDinF01 = Sn'*A1L1*sum(F0TDinF01,2) + mf00TDinF01;
# clear matrices that are not required until the next loop of r
            clear TDin_F F0TDinF* F1TDinF*
# calculate the changes in VALUE ADDED as a weighted average of all TDin-related terms - ADDITIVE
            dv_TDin_F = (v1'-v111TDinF11)/18 + (v110TDinF11-v11011)/18 + (v11100-v111TDinF00)/18 + (v110TDinF00-v11000)/18 + (v11110-v111TDinF10)/36 + (v110TDinF10-v11010)/36 + (v11101-v111TDinF01)/36 + (v110TDinF01-v11001)/36 + (v00111-v001TDinF11)/18 + (v000TDinF11-v00011)/18 + (v00100-v001TDinF00)/18 + (v000TDinF00-v0')/18 + (v00110-v001TDinF10)/36 + (v000TDinF10-v00010)/36 + (v00101-v001TDinF01)/36 + (v000TDinF01-v00001)/36 + (v10111-v101TDinF11)/36 + (v100TDinF11-v10011)/36 + (v10100-v101TDinF00)/36 + (v100TDinF00-v10000)/36 + (v10110-v101TDinF10)/72 + (v100TDinF10-v10010)/72 + (v10101-v101TDinF01)/72 + (v100TDinF01-v10001)/72 + (v01111-v011TDinF11)/36 + (v010TDinF11-v01011)/36 + (v01100-v011TDinF00)/36 + (v010TDinF00-v01000)/36 + (v01110-v011TDinF10)/72 + (v010TDinF10-v01010)/72 + (v01101-v011TDinF01)/72 + (v010TDinF01-v01001)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TDin-related terms - MULTIPLICATIVE
            rv_TDin_F = (((v1'./v111TDinF11) .* (v110TDinF11./v11011) .* (v11100./v111TDinF00) .* (v110TDinF00./v11000)).^(1/18)) .* (((v11110./v111TDinF10) .* (v110TDinF10./v11010) .* (v11101./v111TDinF01) .* (v110TDinF01./v11001)).^(1/36)) .* (((v00111./v001TDinF11) .* (v000TDinF11./v00011) .* (v00100./v001TDinF00) .* (v000TDinF00./v0')).^(1/18)) .* (((v00110./v001TDinF10) .* (v000TDinF10./v00010) .* (v00101./v001TDinF01) .* (v000TDinF01./v00001)).^(1/36)) .* (((v10111./v101TDinF11) .* (v100TDinF11./v10011) .* (v10100./v101TDinF00) .* (v100TDinF00./v10000)).^(1/36)) .* (((v10110./v101TDinF10) .* (v100TDinF10./v10010) .* (v10101./v101TDinF01) .* (v100TDinF01./v10001)).^(1/72)) .* (((v01111./v011TDinF11) .* (v010TDinF11./v01011) .* (v01100./v011TDinF00) .* (v010TDinF00./v01000)).^(1/36)) .* (((v01110./v011TDinF10) .* (v010TDinF10./v01010) .* (v01101./v011TDinF01) .* (v010TDinF01./v01001)).^(1/72));
# store results in structures
            dv.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = dv_TDin_F;
            rv.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = rv_TDin_F;
# calculate the changes in EMPLOYMENT as a weighted average of all TDin-related terms - ADDITIVE
            de_TDin_F = (e1-e111TDinF11)/18 + (e110TDinF11-e11011)/18 + (e11100-e111TDinF00)/18 + (e110TDinF00-e11000)/18 + (e11110-e111TDinF10)/36 + (e110TDinF10-e11010)/36 + (e11101-e111TDinF01)/36 + (e110TDinF01-e11001)/36 + (e00111-e001TDinF11)/18 + (e000TDinF11-e00011)/18 + (e00100-e001TDinF00)/18 + (e000TDinF00-e0)/18 + (e00110-e001TDinF10)/36 + (e000TDinF10-e00010)/36 + (e00101-e001TDinF01)/36 + (e000TDinF01-e00001)/36 + (e10111-e101TDinF11)/36 + (e100TDinF11-e10011)/36 + (e10100-e101TDinF00)/36 + (e100TDinF00-e10000)/36 + (e10110-e101TDinF10)/72 + (e100TDinF10-e10010)/72 + (e10101-e101TDinF01)/72 + (e100TDinF01-e10001)/72 + (e01111-e011TDinF11)/36 + (e010TDinF11-e01011)/36 + (e01100-e011TDinF00)/36 + (e010TDinF00-e01000)/36 + (e01110-e011TDinF10)/72 + (e010TDinF10-e01010)/72 + (e01101-e011TDinF01)/72 + (e010TDinF01-e01001)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TDin-related terms - MULTIPLICATIVE
            re_TDin_F = (((e1./e111TDinF11) .* (e110TDinF11./e11011) .* (e11100./e111TDinF00) .* (e110TDinF00./e11000)).^(1/18)) .* (((e11110./e111TDinF10) .* (e110TDinF10./e11010) .* (e11101./e111TDinF01) .* (e110TDinF01./e11001)).^(1/36)) .* (((e00111./e001TDinF11) .* (e000TDinF11./e00011) .* (e00100./e001TDinF00) .* (e000TDinF00./e0)).^(1/18)) .* (((e00110./e001TDinF10) .* (e000TDinF10./e00010) .* (e00101./e001TDinF01) .* (e000TDinF01./e00001)).^(1/36)) .* (((e10111./e101TDinF11) .* (e100TDinF11./e10011) .* (e10100./e101TDinF00) .* (e100TDinF00./e10000)).^(1/36)) .* (((e10110./e101TDinF10) .* (e100TDinF10./e10010) .* (e10101./e101TDinF01) .* (e100TDinF01./e10001)).^(1/72)) .* (((e01111./e011TDinF11) .* (e010TDinF11./e01011) .* (e01100./e011TDinF00) .* (e010TDinF00./e01000)).^(1/36)) .* (((e01110./e011TDinF10) .* (e010TDinF10./e01010) .* (e01101./e011TDinF01) .* (e010TDinF01./e01001)).^(1/72));
# store results in structures
            de.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = de_TDin_F;
            re.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = re_TDin_F;
# calculate the changes in GDP as a weighted average of all TDin-related terms - ADDITIVE
            dy_TDin_F = (y11111-y111TDinF11)/18 + (y110TDinF11-y11011)/18 + (y11100-y111TDinF00)/18 + (y110TDinF00-y11000)/18 + (y11110-y111TDinF10)/36 + (y110TDinF10-y11010)/36 + (y11101-y111TDinF01)/36 + (y110TDinF01-y11001)/36 + (y00111-y001TDinF11)/18 + (y000TDinF11-y00011)/18 + (y00100-y001TDinF00)/18 + (y000TDinF00-y00000)/18 + (y00110-y001TDinF10)/36 + (y000TDinF10-y00010)/36 + (y00101-y001TDinF01)/36 + (y000TDinF01-y00001)/36 + (y10111-y101TDinF11)/36 + (y100TDinF11-y10011)/36 + (y10100-y101TDinF00)/36 + (y100TDinF00-y10000)/36 + (y10110-y101TDinF10)/72 + (y100TDinF10-y10010)/72 + (y10101-y101TDinF01)/72 + (y100TDinF01-y10001)/72 + (y01111-y011TDinF11)/36 + (y010TDinF11-y01011)/36 + (y01100-y011TDinF00)/36 + (y010TDinF00-y01000)/36 + (y01110-y011TDinF10)/72 + (y010TDinF10-y01010)/72 + (y01101-y011TDinF01)/72 + (y010TDinF01-y01001)/72;
# calculate the changes in GDP as a weighted average of all TDin-related terms - MULTIPLICATIVE
            ry_TDin_F = (((y11111./y111TDinF11) .* (y110TDinF11./y11011) .* (y11100./y111TDinF00) .* (y110TDinF00./y11000)).^(1/18)) .* (((y11110./y111TDinF10) .* (y110TDinF10./y11010) .* (y11101./y111TDinF01) .* (y110TDinF01./y11001)).^(1/36)) .* (((y00111./y001TDinF11) .* (y000TDinF11./y00011) .* (y00100./y001TDinF00) .* (y000TDinF00./y00000)).^(1/18)) .* (((y00110./y001TDinF10) .* (y000TDinF10./y00010) .* (y00101./y001TDinF01) .* (y000TDinF01./y00001)).^(1/36)) .* (((y10111./y101TDinF11) .* (y100TDinF11./y10011) .* (y10100./y101TDinF00) .* (y100TDinF00./y10000)).^(1/36)) .* (((y10110./y101TDinF10) .* (y100TDinF10./y10010) .* (y10101./y101TDinF01) .* (y100TDinF01./y10001)).^(1/72)) .* (((y01111./y011TDinF11) .* (y010TDinF11./y01011) .* (y01100./y011TDinF00) .* (y010TDinF00./y01000)).^(1/36)) .* (((y01110./y011TDinF10) .* (y010TDinF10./y01010) .* (y01101./y011TDinF01) .* (y010TDinF01./y01001)).^(1/72));
# store results in structures
            dy.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = dy_TDin_F;
            ry.TDin.F.(["r",num2str(r),"_s",num2str(s)]) = ry_TDin_F;

# the rest of the (bilateral) component matrices for SDA - TDout
            F1TDoutF11 = (Fcou1-TDout_F).*(Sk'*Find1*diag(f1));
            F0TDoutF11 = (Fcou0+TDout_F).*(Sk'*Find1*diag(f1));
            F1TDoutF00 = (Fcou1-TDout_F).*(Sk'*Find0*diag(f0));
            F0TDoutF00 = (Fcou0+TDout_F).*(Sk'*Find0*diag(f0));
            F1TDoutF10 = (Fcou1-TDout_F).*(Sk'*Find1*diag(f0));
            F0TDoutF10 = (Fcou0+TDout_F).*(Sk'*Find1*diag(f0));
            F1TDoutF01 = (Fcou1-TDout_F).*(Sk'*Find0*diag(f1));
            F0TDoutF01 = (Fcou0+TDout_F).*(Sk'*Find0*diag(f1));
# the rest 32 VALUE ADDED terms for changes in the country of orgin of FINAL demand
# v***TDoutF** = v(Tz*,A*,Fcou*,Find*,f*)
            v111TDoutF11 = Ti1A1L1*sum(F1TDoutF11,2);
            v110TDoutF11 = Ti1A1L1*sum(F0TDoutF11,2);
            v111TDoutF00 = Ti1A1L1*sum(F1TDoutF00,2);
            v110TDoutF00 = Ti1A1L1*sum(F0TDoutF00,2);
            v111TDoutF10 = Ti1A1L1*sum(F1TDoutF10,2);
            v110TDoutF10 = Ti1A1L1*sum(F0TDoutF10,2);
            v111TDoutF01 = Ti1A1L1*sum(F1TDoutF01,2);
            v110TDoutF01 = Ti1A1L1*sum(F0TDoutF01,2);
            v001TDoutF11 = Ti0A0L0*sum(F1TDoutF11,2);
            v000TDoutF11 = Ti0A0L0*sum(F0TDoutF11,2);
            v001TDoutF00 = Ti0A0L0*sum(F1TDoutF00,2);
            v000TDoutF00 = Ti0A0L0*sum(F0TDoutF00,2);
            v001TDoutF10 = Ti0A0L0*sum(F1TDoutF10,2);
            v000TDoutF10 = Ti0A0L0*sum(F0TDoutF10,2);
            v001TDoutF01 = Ti0A0L0*sum(F1TDoutF01,2);
            v000TDoutF01 = Ti0A0L0*sum(F0TDoutF01,2);
            v101TDoutF11 = Ti1A0L0*sum(F1TDoutF11,2);
            v100TDoutF11 = Ti1A0L0*sum(F0TDoutF11,2);
            v101TDoutF00 = Ti1A0L0*sum(F1TDoutF00,2);
            v100TDoutF00 = Ti1A0L0*sum(F0TDoutF00,2);
            v101TDoutF10 = Ti1A0L0*sum(F1TDoutF10,2);
            v100TDoutF10 = Ti1A0L0*sum(F0TDoutF10,2);
            v101TDoutF01 = Ti1A0L0*sum(F1TDoutF01,2);
            v100TDoutF01 = Ti1A0L0*sum(F0TDoutF01,2);
            v011TDoutF11 = Ti0A1L1*sum(F1TDoutF11,2);
            v010TDoutF11 = Ti0A1L1*sum(F0TDoutF11,2);
            v011TDoutF00 = Ti0A1L1*sum(F1TDoutF00,2);
            v010TDoutF00 = Ti0A1L1*sum(F0TDoutF00,2);
            v011TDoutF10 = Ti0A1L1*sum(F1TDoutF10,2);
            v010TDoutF10 = Ti0A1L1*sum(F0TDoutF10,2);
            v011TDoutF01 = Ti0A1L1*sum(F1TDoutF01,2);
            v010TDoutF01 = Ti0A1L1*sum(F0TDoutF01,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of FINAL demand
# e***TDoutF** = e(ec*,A*,Fcou*,Find*,f*)
            e111TDoutF11 = ec1L1*sum(F1TDoutF11,2);
            e110TDoutF11 = ec1L1*sum(F0TDoutF11,2);
            e111TDoutF00 = ec1L1*sum(F1TDoutF00,2);
            e110TDoutF00 = ec1L1*sum(F0TDoutF00,2);
            e111TDoutF10 = ec1L1*sum(F1TDoutF10,2);
            e110TDoutF10 = ec1L1*sum(F0TDoutF10,2);
            e111TDoutF01 = ec1L1*sum(F1TDoutF01,2);
            e110TDoutF01 = ec1L1*sum(F0TDoutF01,2);
            e001TDoutF11 = ec0L0*sum(F1TDoutF11,2);
            e000TDoutF11 = ec0L0*sum(F0TDoutF11,2);
            e001TDoutF00 = ec0L0*sum(F1TDoutF00,2);
            e000TDoutF00 = ec0L0*sum(F0TDoutF00,2);
            e001TDoutF10 = ec0L0*sum(F1TDoutF10,2);
            e000TDoutF10 = ec0L0*sum(F0TDoutF10,2);
            e001TDoutF01 = ec0L0*sum(F1TDoutF01,2);
            e000TDoutF01 = ec0L0*sum(F0TDoutF01,2);
            e101TDoutF11 = ec1L0*sum(F1TDoutF11,2);
            e100TDoutF11 = ec1L0*sum(F0TDoutF11,2);
            e101TDoutF00 = ec1L0*sum(F1TDoutF00,2);
            e100TDoutF00 = ec1L0*sum(F0TDoutF00,2);
            e101TDoutF10 = ec1L0*sum(F1TDoutF10,2);
            e100TDoutF10 = ec1L0*sum(F0TDoutF10,2);
            e101TDoutF01 = ec1L0*sum(F1TDoutF01,2);
            e100TDoutF01 = ec1L0*sum(F0TDoutF01,2);
            e011TDoutF11 = ec0L1*sum(F1TDoutF11,2);
            e010TDoutF11 = ec0L1*sum(F0TDoutF11,2);
            e011TDoutF00 = ec0L1*sum(F1TDoutF00,2);
            e010TDoutF00 = ec0L1*sum(F0TDoutF00,2);
            e011TDoutF10 = ec0L1*sum(F1TDoutF10,2);
            e010TDoutF10 = ec0L1*sum(F0TDoutF10,2);
            e011TDoutF01 = ec0L1*sum(F1TDoutF01,2);
            e010TDoutF01 = ec0L1*sum(F0TDoutF01,2);
# the rest 16 TAXES LESS SUBSIDIES terms for changes in the country of origin of FINAL products
# mf**TDoutF** = mf(Tf*,Fcou*,Find*,f*)
            mf11TDoutF11 = sum(Tf1.*(Sn'*F1TDoutF11),1)';
            mf10TDoutF11 = sum(Tf1.*(Sn'*F0TDoutF11),1)';
            mf11TDoutF00 = sum(Tf1.*(Sn'*F1TDoutF00),1)';
            mf10TDoutF00 = sum(Tf1.*(Sn'*F0TDoutF00),1)';
            mf11TDoutF10 = sum(Tf1.*(Sn'*F1TDoutF10),1)';
            mf10TDoutF10 = sum(Tf1.*(Sn'*F0TDoutF10),1)';
            mf11TDoutF01 = sum(Tf1.*(Sn'*F1TDoutF01),1)';
            mf10TDoutF01 = sum(Tf1.*(Sn'*F0TDoutF01),1)';
            mf01TDoutF11 = sum(Tf0.*(Sn'*F1TDoutF11),1)';
            mf00TDoutF11 = sum(Tf0.*(Sn'*F0TDoutF11),1)';
            mf01TDoutF00 = sum(Tf0.*(Sn'*F1TDoutF00),1)';
            mf00TDoutF00 = sum(Tf0.*(Sn'*F0TDoutF00),1)';
            mf01TDoutF10 = sum(Tf0.*(Sn'*F1TDoutF10),1)';
            mf00TDoutF10 = sum(Tf0.*(Sn'*F0TDoutF10),1)';
            mf01TDoutF01 = sum(Tf0.*(Sn'*F1TDoutF01),1)';
            mf00TDoutF01 = sum(Tf0.*(Sn'*F0TDoutF01),1)';
# the rest 32 GDP terms for changes in the country of origin of FINAL products
# y***TDoutF** = y(Tf*,A,Fcou*,Find*,f*)
            y111TDoutF11 = Sn'*A1L1*sum(F1TDoutF11,2) + mf11TDoutF11;
            y110TDoutF11 = Sn'*A1L1*sum(F0TDoutF11,2) + mf10TDoutF11;
            y111TDoutF00 = Sn'*A1L1*sum(F1TDoutF00,2) + mf11TDoutF00;
            y110TDoutF00 = Sn'*A1L1*sum(F0TDoutF00,2) + mf10TDoutF00;
            y111TDoutF10 = Sn'*A1L1*sum(F1TDoutF10,2) + mf11TDoutF10;
            y110TDoutF10 = Sn'*A1L1*sum(F0TDoutF10,2) + mf10TDoutF10;
            y111TDoutF01 = Sn'*A1L1*sum(F1TDoutF01,2) + mf11TDoutF01;
            y110TDoutF01 = Sn'*A1L1*sum(F0TDoutF01,2) + mf10TDoutF01;
            y001TDoutF11 = Sn'*A0L0*sum(F1TDoutF11,2) + mf01TDoutF11;
            y000TDoutF11 = Sn'*A0L0*sum(F0TDoutF11,2) + mf00TDoutF11;
            y001TDoutF00 = Sn'*A0L0*sum(F1TDoutF00,2) + mf01TDoutF00;
            y000TDoutF00 = Sn'*A0L0*sum(F0TDoutF00,2) + mf00TDoutF00;
            y001TDoutF10 = Sn'*A0L0*sum(F1TDoutF10,2) + mf01TDoutF10;
            y000TDoutF10 = Sn'*A0L0*sum(F0TDoutF10,2) + mf00TDoutF10;
            y001TDoutF01 = Sn'*A0L0*sum(F1TDoutF01,2) + mf01TDoutF01;
            y000TDoutF01 = Sn'*A0L0*sum(F0TDoutF01,2) + mf00TDoutF01;
            y101TDoutF11 = Sn'*A0L0*sum(F1TDoutF11,2) + mf11TDoutF11;
            y100TDoutF11 = Sn'*A0L0*sum(F0TDoutF11,2) + mf10TDoutF11;
            y101TDoutF00 = Sn'*A0L0*sum(F1TDoutF00,2) + mf11TDoutF00;
            y100TDoutF00 = Sn'*A0L0*sum(F0TDoutF00,2) + mf10TDoutF00;
            y101TDoutF10 = Sn'*A0L0*sum(F1TDoutF10,2) + mf11TDoutF10;
            y100TDoutF10 = Sn'*A0L0*sum(F0TDoutF10,2) + mf10TDoutF10;
            y101TDoutF01 = Sn'*A0L0*sum(F1TDoutF01,2) + mf11TDoutF01;
            y100TDoutF01 = Sn'*A0L0*sum(F0TDoutF01,2) + mf10TDoutF01;
            y011TDoutF11 = Sn'*A1L1*sum(F1TDoutF11,2) + mf01TDoutF11;
            y010TDoutF11 = Sn'*A1L1*sum(F0TDoutF11,2) + mf00TDoutF11;
            y011TDoutF00 = Sn'*A1L1*sum(F1TDoutF00,2) + mf01TDoutF00;
            y010TDoutF00 = Sn'*A1L1*sum(F0TDoutF00,2) + mf00TDoutF00;
            y011TDoutF10 = Sn'*A1L1*sum(F1TDoutF10,2) + mf01TDoutF10;
            y010TDoutF10 = Sn'*A1L1*sum(F0TDoutF10,2) + mf00TDoutF10;
            y011TDoutF01 = Sn'*A1L1*sum(F1TDoutF01,2) + mf01TDoutF01;
            y010TDoutF01 = Sn'*A1L1*sum(F0TDoutF01,2) + mf00TDoutF01;
# clear matrices that are not required until the next loop of r
            clear TDout_F F0TDoutF* F1TDoutF*
# calculate the changes in VALUE ADDED as a weighted average of all TDout-related terms - ADDITIVE
            dv_TDout_F = (v1'-v111TDoutF11)/18 + (v110TDoutF11-v11011)/18 + (v11100-v111TDoutF00)/18 + (v110TDoutF00-v11000)/18 + (v11110-v111TDoutF10)/36 + (v110TDoutF10-v11010)/36 + (v11101-v111TDoutF01)/36 + (v110TDoutF01-v11001)/36 + (v00111-v001TDoutF11)/18 + (v000TDoutF11-v00011)/18 + (v00100-v001TDoutF00)/18 + (v000TDoutF00-v0')/18 + (v00110-v001TDoutF10)/36 + (v000TDoutF10-v00010)/36 + (v00101-v001TDoutF01)/36 + (v000TDoutF01-v00001)/36 + (v10111-v101TDoutF11)/36 + (v100TDoutF11-v10011)/36 + (v10100-v101TDoutF00)/36 + (v100TDoutF00-v10000)/36 + (v10110-v101TDoutF10)/72 + (v100TDoutF10-v10010)/72 + (v10101-v101TDoutF01)/72 + (v100TDoutF01-v10001)/72 + (v01111-v011TDoutF11)/36 + (v010TDoutF11-v01011)/36 + (v01100-v011TDoutF00)/36 + (v010TDoutF00-v01000)/36 + (v01110-v011TDoutF10)/72 + (v010TDoutF10-v01010)/72 + (v01101-v011TDoutF01)/72 + (v010TDoutF01-v01001)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TDout-related terms - MULTIPLICATIVE
            rv_TDout_F = (((v1'./v111TDoutF11) .* (v110TDoutF11./v11011) .* (v11100./v111TDoutF00) .* (v110TDoutF00./v11000)).^(1/18)) .* (((v11110./v111TDoutF10) .* (v110TDoutF10./v11010) .* (v11101./v111TDoutF01) .* (v110TDoutF01./v11001)).^(1/36)) .* (((v00111./v001TDoutF11) .* (v000TDoutF11./v00011) .* (v00100./v001TDoutF00) .* (v000TDoutF00./v0')).^(1/18)) .* (((v00110./v001TDoutF10) .* (v000TDoutF10./v00010) .* (v00101./v001TDoutF01) .* (v000TDoutF01./v00001)).^(1/36)) .* (((v10111./v101TDoutF11) .* (v100TDoutF11./v10011) .* (v10100./v101TDoutF00) .* (v100TDoutF00./v10000)).^(1/36)) .* (((v10110./v101TDoutF10) .* (v100TDoutF10./v10010) .* (v10101./v101TDoutF01) .* (v100TDoutF01./v10001)).^(1/72)) .* (((v01111./v011TDoutF11) .* (v010TDoutF11./v01011) .* (v01100./v011TDoutF00) .* (v010TDoutF00./v01000)).^(1/36)) .* (((v01110./v011TDoutF10) .* (v010TDoutF10./v01010) .* (v01101./v011TDoutF01) .* (v010TDoutF01./v01001)).^(1/72));
# store results in structures
            dv.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = dv_TDout_F;
            rv.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = rv_TDout_F;
# calculate the changes in EMPLOYMENT as a weighted average of all TDout-related terms - ADDITIVE
            de_TDout_F = (e1-e111TDoutF11)/18 + (e110TDoutF11-e11011)/18 + (e11100-e111TDoutF00)/18 + (e110TDoutF00-e11000)/18 + (e11110-e111TDoutF10)/36 + (e110TDoutF10-e11010)/36 + (e11101-e111TDoutF01)/36 + (e110TDoutF01-e11001)/36 + (e00111-e001TDoutF11)/18 + (e000TDoutF11-e00011)/18 + (e00100-e001TDoutF00)/18 + (e000TDoutF00-e0)/18 + (e00110-e001TDoutF10)/36 + (e000TDoutF10-e00010)/36 + (e00101-e001TDoutF01)/36 + (e000TDoutF01-e00001)/36 + (e10111-e101TDoutF11)/36 + (e100TDoutF11-e10011)/36 + (e10100-e101TDoutF00)/36 + (e100TDoutF00-e10000)/36 + (e10110-e101TDoutF10)/72 + (e100TDoutF10-e10010)/72 + (e10101-e101TDoutF01)/72 + (e100TDoutF01-e10001)/72 + (e01111-e011TDoutF11)/36 + (e010TDoutF11-e01011)/36 + (e01100-e011TDoutF00)/36 + (e010TDoutF00-e01000)/36 + (e01110-e011TDoutF10)/72 + (e010TDoutF10-e01010)/72 + (e01101-e011TDoutF01)/72 + (e010TDoutF01-e01001)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TDout-related terms - MULTIPLICATIVE
            re_TDout_F = (((e1./e111TDoutF11) .* (e110TDoutF11./e11011) .* (e11100./e111TDoutF00) .* (e110TDoutF00./e11000)).^(1/18)) .* (((e11110./e111TDoutF10) .* (e110TDoutF10./e11010) .* (e11101./e111TDoutF01) .* (e110TDoutF01./e11001)).^(1/36)) .* (((e00111./e001TDoutF11) .* (e000TDoutF11./e00011) .* (e00100./e001TDoutF00) .* (e000TDoutF00./e0)).^(1/18)) .* (((e00110./e001TDoutF10) .* (e000TDoutF10./e00010) .* (e00101./e001TDoutF01) .* (e000TDoutF01./e00001)).^(1/36)) .* (((e10111./e101TDoutF11) .* (e100TDoutF11./e10011) .* (e10100./e101TDoutF00) .* (e100TDoutF00./e10000)).^(1/36)) .* (((e10110./e101TDoutF10) .* (e100TDoutF10./e10010) .* (e10101./e101TDoutF01) .* (e100TDoutF01./e10001)).^(1/72)) .* (((e01111./e011TDoutF11) .* (e010TDoutF11./e01011) .* (e01100./e011TDoutF00) .* (e010TDoutF00./e01000)).^(1/36)) .* (((e01110./e011TDoutF10) .* (e010TDoutF10./e01010) .* (e01101./e011TDoutF01) .* (e010TDoutF01./e01001)).^(1/72));
# store results in structures
            de.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = de_TDout_F;
            re.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = re_TDout_F;
# calculate the changes in GDP as a weighted average of all TDout-related terms - ADDITIVE
            dy_TDout_F = (y11111-y111TDoutF11)/18 + (y110TDoutF11-y11011)/18 + (y11100-y111TDoutF00)/18 + (y110TDoutF00-y11000)/18 + (y11110-y111TDoutF10)/36 + (y110TDoutF10-y11010)/36 + (y11101-y111TDoutF01)/36 + (y110TDoutF01-y11001)/36 + (y00111-y001TDoutF11)/18 + (y000TDoutF11-y00011)/18 + (y00100-y001TDoutF00)/18 + (y000TDoutF00-y00000)/18 + (y00110-y001TDoutF10)/36 + (y000TDoutF10-y00010)/36 + (y00101-y001TDoutF01)/36 + (y000TDoutF01-y00001)/36 + (y10111-y101TDoutF11)/36 + (y100TDoutF11-y10011)/36 + (y10100-y101TDoutF00)/36 + (y100TDoutF00-y10000)/36 + (y10110-y101TDoutF10)/72 + (y100TDoutF10-y10010)/72 + (y10101-y101TDoutF01)/72 + (y100TDoutF01-y10001)/72 + (y01111-y011TDoutF11)/36 + (y010TDoutF11-y01011)/36 + (y01100-y011TDoutF00)/36 + (y010TDoutF00-y01000)/36 + (y01110-y011TDoutF10)/72 + (y010TDoutF10-y01010)/72 + (y01101-y011TDoutF01)/72 + (y010TDoutF01-y01001)/72;
# calculate the changes in GDP as a weighted average of all TDout-related terms - MULTIPLICATIVE
            ry_TDout_F = (((y11111./y111TDoutF11) .* (y110TDoutF11./y11011) .* (y11100./y111TDoutF00) .* (y110TDoutF00./y11000)).^(1/18)) .* (((y11110./y111TDoutF10) .* (y110TDoutF10./y11010) .* (y11101./y111TDoutF01) .* (y110TDoutF01./y11001)).^(1/36)) .* (((y00111./y001TDoutF11) .* (y000TDoutF11./y00011) .* (y00100./y001TDoutF00) .* (y000TDoutF00./y00000)).^(1/18)) .* (((y00110./y001TDoutF10) .* (y000TDoutF10./y00010) .* (y00101./y001TDoutF01) .* (y000TDoutF01./y00001)).^(1/36)) .* (((y10111./y101TDoutF11) .* (y100TDoutF11./y10011) .* (y10100./y101TDoutF00) .* (y100TDoutF00./y10000)).^(1/36)) .* (((y10110./y101TDoutF10) .* (y100TDoutF10./y10010) .* (y10101./y101TDoutF01) .* (y100TDoutF01./y10001)).^(1/72)) .* (((y01111./y011TDoutF11) .* (y010TDoutF11./y01011) .* (y01100./y011TDoutF00) .* (y010TDoutF00./y01000)).^(1/36)) .* (((y01110./y011TDoutF10) .* (y010TDoutF10./y01010) .* (y01101./y011TDoutF01) .* (y010TDoutF01./y01001)).^(1/72));
# store results in structures
            dy.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = dy_TDout_F;
            ry.TDout.F.(["r",num2str(r),"_s",num2str(s)]) = ry_TDout_F;

# the rest of the (bilateral) component matrices for SDA - TR
            F1TRF11 = (Fcou1-TR_F).*(Sk'*Find1*diag(f1));
            F0TRF11 = (Fcou0+TR_F).*(Sk'*Find1*diag(f1));
            F1TRF00 = (Fcou1-TR_F).*(Sk'*Find0*diag(f0));
            F0TRF00 = (Fcou0+TR_F).*(Sk'*Find0*diag(f0));
            F1TRF10 = (Fcou1-TR_F).*(Sk'*Find1*diag(f0));
            F0TRF10 = (Fcou0+TR_F).*(Sk'*Find1*diag(f0));
            F1TRF01 = (Fcou1-TR_F).*(Sk'*Find0*diag(f1));
            F0TRF01 = (Fcou0+TR_F).*(Sk'*Find0*diag(f1));
# the rest 32 VALUE ADDED terms for changes in the country of orgin of FINAL demand
# v***TRF** = v(Tz*,A*,Fcou*,Find*,f*)
            v111TRF11 = Ti1A1L1*sum(F1TRF11,2);
            v110TRF11 = Ti1A1L1*sum(F0TRF11,2);
            v111TRF00 = Ti1A1L1*sum(F1TRF00,2);
            v110TRF00 = Ti1A1L1*sum(F0TRF00,2);
            v111TRF10 = Ti1A1L1*sum(F1TRF10,2);
            v110TRF10 = Ti1A1L1*sum(F0TRF10,2);
            v111TRF01 = Ti1A1L1*sum(F1TRF01,2);
            v110TRF01 = Ti1A1L1*sum(F0TRF01,2);
            v001TRF11 = Ti0A0L0*sum(F1TRF11,2);
            v000TRF11 = Ti0A0L0*sum(F0TRF11,2);
            v001TRF00 = Ti0A0L0*sum(F1TRF00,2);
            v000TRF00 = Ti0A0L0*sum(F0TRF00,2);
            v001TRF10 = Ti0A0L0*sum(F1TRF10,2);
            v000TRF10 = Ti0A0L0*sum(F0TRF10,2);
            v001TRF01 = Ti0A0L0*sum(F1TRF01,2);
            v000TRF01 = Ti0A0L0*sum(F0TRF01,2);
            v101TRF11 = Ti1A0L0*sum(F1TRF11,2);
            v100TRF11 = Ti1A0L0*sum(F0TRF11,2);
            v101TRF00 = Ti1A0L0*sum(F1TRF00,2);
            v100TRF00 = Ti1A0L0*sum(F0TRF00,2);
            v101TRF10 = Ti1A0L0*sum(F1TRF10,2);
            v100TRF10 = Ti1A0L0*sum(F0TRF10,2);
            v101TRF01 = Ti1A0L0*sum(F1TRF01,2);
            v100TRF01 = Ti1A0L0*sum(F0TRF01,2);
            v011TRF11 = Ti0A1L1*sum(F1TRF11,2);
            v010TRF11 = Ti0A1L1*sum(F0TRF11,2);
            v011TRF00 = Ti0A1L1*sum(F1TRF00,2);
            v010TRF00 = Ti0A1L1*sum(F0TRF00,2);
            v011TRF10 = Ti0A1L1*sum(F1TRF10,2);
            v010TRF10 = Ti0A1L1*sum(F0TRF10,2);
            v011TRF01 = Ti0A1L1*sum(F1TRF01,2);
            v010TRF01 = Ti0A1L1*sum(F0TRF01,2);
# the rest 32 EMPLOYMENT terms for changes in the country of orgin of FINAL demand
# e***TRF** = e(ec*,A*,Fcou*,Find*,f*)
            e111TRF11 = ec1L1*sum(F1TRF11,2);
            e110TRF11 = ec1L1*sum(F0TRF11,2);
            e111TRF00 = ec1L1*sum(F1TRF00,2);
            e110TRF00 = ec1L1*sum(F0TRF00,2);
            e111TRF10 = ec1L1*sum(F1TRF10,2);
            e110TRF10 = ec1L1*sum(F0TRF10,2);
            e111TRF01 = ec1L1*sum(F1TRF01,2);
            e110TRF01 = ec1L1*sum(F0TRF01,2);
            e001TRF11 = ec0L0*sum(F1TRF11,2);
            e000TRF11 = ec0L0*sum(F0TRF11,2);
            e001TRF00 = ec0L0*sum(F1TRF00,2);
            e000TRF00 = ec0L0*sum(F0TRF00,2);
            e001TRF10 = ec0L0*sum(F1TRF10,2);
            e000TRF10 = ec0L0*sum(F0TRF10,2);
            e001TRF01 = ec0L0*sum(F1TRF01,2);
            e000TRF01 = ec0L0*sum(F0TRF01,2);
            e101TRF11 = ec1L0*sum(F1TRF11,2);
            e100TRF11 = ec1L0*sum(F0TRF11,2);
            e101TRF00 = ec1L0*sum(F1TRF00,2);
            e100TRF00 = ec1L0*sum(F0TRF00,2);
            e101TRF10 = ec1L0*sum(F1TRF10,2);
            e100TRF10 = ec1L0*sum(F0TRF10,2);
            e101TRF01 = ec1L0*sum(F1TRF01,2);
            e100TRF01 = ec1L0*sum(F0TRF01,2);
            e011TRF11 = ec0L1*sum(F1TRF11,2);
            e010TRF11 = ec0L1*sum(F0TRF11,2);
            e011TRF00 = ec0L1*sum(F1TRF00,2);
            e010TRF00 = ec0L1*sum(F0TRF00,2);
            e011TRF10 = ec0L1*sum(F1TRF10,2);
            e010TRF10 = ec0L1*sum(F0TRF10,2);
            e011TRF01 = ec0L1*sum(F1TRF01,2);
            e010TRF01 = ec0L1*sum(F0TRF01,2);
# the rest 16 TAXES LESS SUBSIDIES terms for changes in the country of origin of FINAL products
# mf**TRF** = mf(Tf*,Fcou*,Find*,f*)
            mf11TRF11 = sum(Tf1.*(Sn'*F1TRF11),1)';
            mf10TRF11 = sum(Tf1.*(Sn'*F0TRF11),1)';
            mf11TRF00 = sum(Tf1.*(Sn'*F1TRF00),1)';
            mf10TRF00 = sum(Tf1.*(Sn'*F0TRF00),1)';
            mf11TRF10 = sum(Tf1.*(Sn'*F1TRF10),1)';
            mf10TRF10 = sum(Tf1.*(Sn'*F0TRF10),1)';
            mf11TRF01 = sum(Tf1.*(Sn'*F1TRF01),1)';
            mf10TRF01 = sum(Tf1.*(Sn'*F0TRF01),1)';
            mf01TRF11 = sum(Tf0.*(Sn'*F1TRF11),1)';
            mf00TRF11 = sum(Tf0.*(Sn'*F0TRF11),1)';
            mf01TRF00 = sum(Tf0.*(Sn'*F1TRF00),1)';
            mf00TRF00 = sum(Tf0.*(Sn'*F0TRF00),1)';
            mf01TRF10 = sum(Tf0.*(Sn'*F1TRF10),1)';
            mf00TRF10 = sum(Tf0.*(Sn'*F0TRF10),1)';
            mf01TRF01 = sum(Tf0.*(Sn'*F1TRF01),1)';
            mf00TRF01 = sum(Tf0.*(Sn'*F0TRF01),1)';
# the rest 32 GDP terms for changes in the country of origin of FINAL products
# y***TRF** = y(Tf*,A,Fcou*,Find*,f*)
            y111TRF11 = Sn'*A1L1*sum(F1TRF11,2) + mf11TRF11;
            y110TRF11 = Sn'*A1L1*sum(F0TRF11,2) + mf10TRF11;
            y111TRF00 = Sn'*A1L1*sum(F1TRF00,2) + mf11TRF00;
            y110TRF00 = Sn'*A1L1*sum(F0TRF00,2) + mf10TRF00;
            y111TRF10 = Sn'*A1L1*sum(F1TRF10,2) + mf11TRF10;
            y110TRF10 = Sn'*A1L1*sum(F0TRF10,2) + mf10TRF10;
            y111TRF01 = Sn'*A1L1*sum(F1TRF01,2) + mf11TRF01;
            y110TRF01 = Sn'*A1L1*sum(F0TRF01,2) + mf10TRF01;
            y001TRF11 = Sn'*A0L0*sum(F1TRF11,2) + mf01TRF11;
            y000TRF11 = Sn'*A0L0*sum(F0TRF11,2) + mf00TRF11;
            y001TRF00 = Sn'*A0L0*sum(F1TRF00,2) + mf01TRF00;
            y000TRF00 = Sn'*A0L0*sum(F0TRF00,2) + mf00TRF00;
            y001TRF10 = Sn'*A0L0*sum(F1TRF10,2) + mf01TRF10;
            y000TRF10 = Sn'*A0L0*sum(F0TRF10,2) + mf00TRF10;
            y001TRF01 = Sn'*A0L0*sum(F1TRF01,2) + mf01TRF01;
            y000TRF01 = Sn'*A0L0*sum(F0TRF01,2) + mf00TRF01;
            y101TRF11 = Sn'*A0L0*sum(F1TRF11,2) + mf11TRF11;
            y100TRF11 = Sn'*A0L0*sum(F0TRF11,2) + mf10TRF11;
            y101TRF00 = Sn'*A0L0*sum(F1TRF00,2) + mf11TRF00;
            y100TRF00 = Sn'*A0L0*sum(F0TRF00,2) + mf10TRF00;
            y101TRF10 = Sn'*A0L0*sum(F1TRF10,2) + mf11TRF10;
            y100TRF10 = Sn'*A0L0*sum(F0TRF10,2) + mf10TRF10;
            y101TRF01 = Sn'*A0L0*sum(F1TRF01,2) + mf11TRF01;
            y100TRF01 = Sn'*A0L0*sum(F0TRF01,2) + mf10TRF01;
            y011TRF11 = Sn'*A1L1*sum(F1TRF11,2) + mf01TRF11;
            y010TRF11 = Sn'*A1L1*sum(F0TRF11,2) + mf00TRF11;
            y011TRF00 = Sn'*A1L1*sum(F1TRF00,2) + mf01TRF00;
            y010TRF00 = Sn'*A1L1*sum(F0TRF00,2) + mf00TRF00;
            y011TRF10 = Sn'*A1L1*sum(F1TRF10,2) + mf01TRF10;
            y010TRF10 = Sn'*A1L1*sum(F0TRF10,2) + mf00TRF10;
            y011TRF01 = Sn'*A1L1*sum(F1TRF01,2) + mf01TRF01;
            y010TRF01 = Sn'*A1L1*sum(F0TRF01,2) + mf00TRF01;
# clear matrices that are not required until the next loop of r
            clear TR_F F0TRF* F1TRF*
# calculate the changes in VALUE ADDED as a weighted average of all TR-related terms - ADDITIVE
            dv_TR_F = (v1'-v111TRF11)/18 + (v110TRF11-v11011)/18 + (v11100-v111TRF00)/18 + (v110TRF00-v11000)/18 + (v11110-v111TRF10)/36 + (v110TRF10-v11010)/36 + (v11101-v111TRF01)/36 + (v110TRF01-v11001)/36 + (v00111-v001TRF11)/18 + (v000TRF11-v00011)/18 + (v00100-v001TRF00)/18 + (v000TRF00-v0')/18 + (v00110-v001TRF10)/36 + (v000TRF10-v00010)/36 + (v00101-v001TRF01)/36 + (v000TRF01-v00001)/36 + (v10111-v101TRF11)/36 + (v100TRF11-v10011)/36 + (v10100-v101TRF00)/36 + (v100TRF00-v10000)/36 + (v10110-v101TRF10)/72 + (v100TRF10-v10010)/72 + (v10101-v101TRF01)/72 + (v100TRF01-v10001)/72 + (v01111-v011TRF11)/36 + (v010TRF11-v01011)/36 + (v01100-v011TRF00)/36 + (v010TRF00-v01000)/36 + (v01110-v011TRF10)/72 + (v010TRF10-v01010)/72 + (v01101-v011TRF01)/72 + (v010TRF01-v01001)/72;
# calculate the changes in VALUE ADDED as a weighted average of all TR-related terms - MULTIPLICATIVE
            rv_TR_F = (((v1'./v111TRF11) .* (v110TRF11./v11011) .* (v11100./v111TRF00) .* (v110TRF00./v11000)).^(1/18)) .* (((v11110./v111TRF10) .* (v110TRF10./v11010) .* (v11101./v111TRF01) .* (v110TRF01./v11001)).^(1/36)) .* (((v00111./v001TRF11) .* (v000TRF11./v00011) .* (v00100./v001TRF00) .* (v000TRF00./v0')).^(1/18)) .* (((v00110./v001TRF10) .* (v000TRF10./v00010) .* (v00101./v001TRF01) .* (v000TRF01./v00001)).^(1/36)) .* (((v10111./v101TRF11) .* (v100TRF11./v10011) .* (v10100./v101TRF00) .* (v100TRF00./v10000)).^(1/36)) .* (((v10110./v101TRF10) .* (v100TRF10./v10010) .* (v10101./v101TRF01) .* (v100TRF01./v10001)).^(1/72)) .* (((v01111./v011TRF11) .* (v010TRF11./v01011) .* (v01100./v011TRF00) .* (v010TRF00./v01000)).^(1/36)) .* (((v01110./v011TRF10) .* (v010TRF10./v01010) .* (v01101./v011TRF01) .* (v010TRF01./v01001)).^(1/72));
# store results in structures
            dv.TR.F.(["r",num2str(r),"_s",num2str(s)]) = dv_TR_F;
            rv.TR.F.(["r",num2str(r),"_s",num2str(s)]) = rv_TR_F;
# calculate the changes in EMPLOYMENT as a weighted average of all TR-related terms - ADDITIVE
            de_TR_F = (e1-e111TRF11)/18 + (e110TRF11-e11011)/18 + (e11100-e111TRF00)/18 + (e110TRF00-e11000)/18 + (e11110-e111TRF10)/36 + (e110TRF10-e11010)/36 + (e11101-e111TRF01)/36 + (e110TRF01-e11001)/36 + (e00111-e001TRF11)/18 + (e000TRF11-e00011)/18 + (e00100-e001TRF00)/18 + (e000TRF00-e0)/18 + (e00110-e001TRF10)/36 + (e000TRF10-e00010)/36 + (e00101-e001TRF01)/36 + (e000TRF01-e00001)/36 + (e10111-e101TRF11)/36 + (e100TRF11-e10011)/36 + (e10100-e101TRF00)/36 + (e100TRF00-e10000)/36 + (e10110-e101TRF10)/72 + (e100TRF10-e10010)/72 + (e10101-e101TRF01)/72 + (e100TRF01-e10001)/72 + (e01111-e011TRF11)/36 + (e010TRF11-e01011)/36 + (e01100-e011TRF00)/36 + (e010TRF00-e01000)/36 + (e01110-e011TRF10)/72 + (e010TRF10-e01010)/72 + (e01101-e011TRF01)/72 + (e010TRF01-e01001)/72;
# calculate the changes in EMPLOYMENT as a weighted average of all TR-related terms - MULTIPLICATIVE
            re_TR_F = (((e1./e111TRF11) .* (e110TRF11./e11011) .* (e11100./e111TRF00) .* (e110TRF00./e11000)).^(1/18)) .* (((e11110./e111TRF10) .* (e110TRF10./e11010) .* (e11101./e111TRF01) .* (e110TRF01./e11001)).^(1/36)) .* (((e00111./e001TRF11) .* (e000TRF11./e00011) .* (e00100./e001TRF00) .* (e000TRF00./e0)).^(1/18)) .* (((e00110./e001TRF10) .* (e000TRF10./e00010) .* (e00101./e001TRF01) .* (e000TRF01./e00001)).^(1/36)) .* (((e10111./e101TRF11) .* (e100TRF11./e10011) .* (e10100./e101TRF00) .* (e100TRF00./e10000)).^(1/36)) .* (((e10110./e101TRF10) .* (e100TRF10./e10010) .* (e10101./e101TRF01) .* (e100TRF01./e10001)).^(1/72)) .* (((e01111./e011TRF11) .* (e010TRF11./e01011) .* (e01100./e011TRF00) .* (e010TRF00./e01000)).^(1/36)) .* (((e01110./e011TRF10) .* (e010TRF10./e01010) .* (e01101./e011TRF01) .* (e010TRF01./e01001)).^(1/72));
# store results in structures
            de.TR.F.(["r",num2str(r),"_s",num2str(s)]) = de_TR_F;
            re.TR.F.(["r",num2str(r),"_s",num2str(s)]) = re_TR_F;
# calculate the changes in GDP as a weighted average of all TR-related terms - ADDITIVE
            dy_TR_F = (y11111-y111TRF11)/18 + (y110TRF11-y11011)/18 + (y11100-y111TRF00)/18 + (y110TRF00-y11000)/18 + (y11110-y111TRF10)/36 + (y110TRF10-y11010)/36 + (y11101-y111TRF01)/36 + (y110TRF01-y11001)/36 + (y00111-y001TRF11)/18 + (y000TRF11-y00011)/18 + (y00100-y001TRF00)/18 + (y000TRF00-y00000)/18 + (y00110-y001TRF10)/36 + (y000TRF10-y00010)/36 + (y00101-y001TRF01)/36 + (y000TRF01-y00001)/36 + (y10111-y101TRF11)/36 + (y100TRF11-y10011)/36 + (y10100-y101TRF00)/36 + (y100TRF00-y10000)/36 + (y10110-y101TRF10)/72 + (y100TRF10-y10010)/72 + (y10101-y101TRF01)/72 + (y100TRF01-y10001)/72 + (y01111-y011TRF11)/36 + (y010TRF11-y01011)/36 + (y01100-y011TRF00)/36 + (y010TRF00-y01000)/36 + (y01110-y011TRF10)/72 + (y010TRF10-y01010)/72 + (y01101-y011TRF01)/72 + (y010TRF01-y01001)/72;
# calculate the changes in GDP as a weighted average of all TR-related terms - MULTIPLICATIVE
            ry_TR_F = (((y11111./y111TRF11) .* (y110TRF11./y11011) .* (y11100./y111TRF00) .* (y110TRF00./y11000)).^(1/18)) .* (((y11110./y111TRF10) .* (y110TRF10./y11010) .* (y11101./y111TRF01) .* (y110TRF01./y11001)).^(1/36)) .* (((y00111./y001TRF11) .* (y000TRF11./y00011) .* (y00100./y001TRF00) .* (y000TRF00./y00000)).^(1/18)) .* (((y00110./y001TRF10) .* (y000TRF10./y00010) .* (y00101./y001TRF01) .* (y000TRF01./y00001)).^(1/36)) .* (((y10111./y101TRF11) .* (y100TRF11./y10011) .* (y10100./y101TRF00) .* (y100TRF00./y10000)).^(1/36)) .* (((y10110./y101TRF10) .* (y100TRF10./y10010) .* (y10101./y101TRF01) .* (y100TRF01./y10001)).^(1/72)) .* (((y01111./y011TRF11) .* (y010TRF11./y01011) .* (y01100./y011TRF00) .* (y010TRF00./y01000)).^(1/36)) .* (((y01110./y011TRF10) .* (y010TRF10./y01010) .* (y01101./y011TRF01) .* (y010TRF01./y01001)).^(1/72));
# store results in structures
            dy.TR.F.(["r",num2str(r),"_s",num2str(s)]) = dy_TR_F;
            ry.TR.F.(["r",num2str(r),"_s",num2str(s)]) = ry_TR_F;

# end the if statement for (r != s)
        endif
# end loop for r
    endfor
# clear matrices that are not required until the next loop of s
    clear TCD_F*
# end loop for s
endfor

# clear unnecessary variables to save space
clear A* L* v00* v11* v01* v10* e00* e11* e01* e10* y001* y0001* y00001 y110* y1110* y11110 y01* y10* mf00* mf11* mf01* mf10*

# print a message to monitor progress (optional)
printf ("All SDA calculations have been completed. Now combining results for intermediate and final demand before writing the results into CSV files.\n")

# combine the results and write to a CSV file
# define file name to extract the selected data (only for the r-s country pair)
filename = ["/[insert the path to the output CSV files]/KORUS_bp_" y0 "-" y1 ".csv"];
# write the first line in the file, so that the variables in loops below are appended to it
varname = ["KORUS_bp_" y0 "-" y1];
dlmwrite (filename, varname, ",")

for s = [19,36]
    for r = [36,19]
# ignore a case where r = s
        if (r != s)
# combine the results for intermediate and final demand - ADDITIVE
        dv.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = dv.TC.A.(["r",num2str(r),"_s",num2str(s)]) + dv.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        dv.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = dv.TDin.A.(["r",num2str(r),"_s",num2str(s)]) + dv.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        dv.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = dv.TDout.A.(["r",num2str(r),"_s",num2str(s)]) + dv.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        dv.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = dv.TR.A.(["r",num2str(r),"_s",num2str(s)]) + dv.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        de.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = de.TC.A.(["r",num2str(r),"_s",num2str(s)]) + de.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        de.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = de.TDin.A.(["r",num2str(r),"_s",num2str(s)]) + de.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        de.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = de.TDout.A.(["r",num2str(r),"_s",num2str(s)]) + de.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        de.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = de.TR.A.(["r",num2str(r),"_s",num2str(s)]) + de.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        dy.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = dy.TC.A.(["r",num2str(r),"_s",num2str(s)]) + dy.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        dy.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = dy.TDin.A.(["r",num2str(r),"_s",num2str(s)]) + dy.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        dy.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = dy.TDout.A.(["r",num2str(r),"_s",num2str(s)]) + dy.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        dy.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = dy.TR.A.(["r",num2str(r),"_s",num2str(s)]) + dy.TR.F.(["r",num2str(r),"_s",num2str(s)]);
# combine the results for intermediate and final demand - MULTIPICATIVE
        rv.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = rv.TC.A.(["r",num2str(r),"_s",num2str(s)]) .* rv.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        rv.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = rv.TDin.A.(["r",num2str(r),"_s",num2str(s)]) .* rv.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        rv.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = rv.TDout.A.(["r",num2str(r),"_s",num2str(s)]) .* rv.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        rv.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = rv.TR.A.(["r",num2str(r),"_s",num2str(s)]) .* rv.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        re.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = re.TC.A.(["r",num2str(r),"_s",num2str(s)]) .* re.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        re.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = re.TDin.A.(["r",num2str(r),"_s",num2str(s)]) .* re.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        re.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = re.TDout.A.(["r",num2str(r),"_s",num2str(s)]) .* re.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        re.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = re.TR.A.(["r",num2str(r),"_s",num2str(s)]) .* re.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        ry.TC.AF.(["r",num2str(r),"_s",num2str(s)]) = ry.TC.A.(["r",num2str(r),"_s",num2str(s)]) .* ry.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        ry.TDin.AF.(["r",num2str(r),"_s",num2str(s)]) = ry.TDin.A.(["r",num2str(r),"_s",num2str(s)]) .* ry.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        ry.TDout.AF.(["r",num2str(r),"_s",num2str(s)]) = ry.TDout.A.(["r",num2str(r),"_s",num2str(s)]) .* ry.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        ry.TR.AF.(["r",num2str(r),"_s",num2str(s)]) = ry.TR.A.(["r",num2str(r),"_s",num2str(s)]) .* ry.TR.F.(["r",num2str(r),"_s",num2str(s)]);
# define variables and write to the CSV file
# each variable [dv*, de*, rv*, re*] is reshaped and transposed into a 65 country x 36 industry matrix
# each variable [dy*, ry*] is kept as a 65 country x 1 column vector
# combined INTERMEDIATE and FINAL demand
# dv.**.AF
        varname = ["dvTCAF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTRAF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# de.**.AF
        varname = ["deTCAF_" "r" num2str(r) "s" num2str(s)];
        var = de.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = de.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = de.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTRAF_" "r" num2str(r) "s" num2str(s)];
        var = de.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# dy.**.AF
        varname = ["dyTCAF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTRAF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")
# rv.**.AF
        varname = ["rvTCAF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTRAF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# re.**.AF
        varname = ["reTCAF_" "r" num2str(r) "s" num2str(s)];
        var = re.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = re.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = re.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTRAF_" "r" num2str(r) "s" num2str(s)];
        var = re.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# ry.**.AF
        varname = ["ryTCAF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TC.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDinAF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDin.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDoutAF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDout.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTRAF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TR.AF.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")
# INTERMEDIATE demand
# dv.**.A
        varname = ["dvTCA_" "r" num2str(r) "s" num2str(s)];
        var = dv.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDinA_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTRA_" "r" num2str(r) "s" num2str(s)];
        var = dv.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# de.**.A
        varname = ["deTCA_" "r" num2str(r) "s" num2str(s)];
        var = de.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDinA_" "r" num2str(r) "s" num2str(s)];
        var = de.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = de.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTRA_" "r" num2str(r) "s" num2str(s)];
        var = de.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# dy.**.A
        varname = ["dyTCA_" "r" num2str(r) "s" num2str(s)];
        var = dy.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDinA_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTRA_" "r" num2str(r) "s" num2str(s)];
        var = dy.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")
# rv.**.A
        varname = ["rvTCA_" "r" num2str(r) "s" num2str(s)];
        var = rv.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDinA_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTRA_" "r" num2str(r) "s" num2str(s)];
        var = rv.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# re.**.A
        varname = ["reTCA_" "r" num2str(r) "s" num2str(s)];
        var = re.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDinA_" "r" num2str(r) "s" num2str(s)];
        var = re.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = re.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTRA_" "r" num2str(r) "s" num2str(s)];
        var = re.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# ry.**.A
        varname = ["ryTCA_" "r" num2str(r) "s" num2str(s)];
        var = ry.TC.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDinA_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDin.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDoutA_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDout.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTRA_" "r" num2str(r) "s" num2str(s)];
        var = ry.TR.A.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")
# FINAL demand
# dv.**.F
        varname = ["dvTCF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDinF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["dvTRF_" "r" num2str(r) "s" num2str(s)];
        var = dv.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# de.**.F
        varname = ["deTCF_" "r" num2str(r) "s" num2str(s)];
        var = de.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDinF_" "r" num2str(r) "s" num2str(s)];
        var = de.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = de.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["deTRF_" "r" num2str(r) "s" num2str(s)];
        var = de.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# dy.**.F
        varname = ["dyTCF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDinF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["dyTRF_" "r" num2str(r) "s" num2str(s)];
        var = dy.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")
# rv.**.F
        varname = ["rvTCF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDinF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["rvTRF_" "r" num2str(r) "s" num2str(s)];
        var = rv.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# re.**.F
        varname = ["reTCF_" "r" num2str(r) "s" num2str(s)];
        var = re.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDinF_" "r" num2str(r) "s" num2str(s)];
        var = re.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = re.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

        varname = ["reTRF_" "r" num2str(r) "s" num2str(s)];
        var = re.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        var = reshape(var,36,65);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")
# ry.**.F
        varname = ["ryTCF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TC.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDinF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDin.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTDoutF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TDout.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

        varname = ["ryTRF_" "r" num2str(r) "s" num2str(s)];
        var = ry.TR.F.(["r",num2str(r),"_s",num2str(s)]);
        dlmwrite (filename, varname, ",", "-append")
        dlmwrite (filename, [var(s);var(r)], ",", "-append")

# end the if statement for (r != s)
        endif
# end loops for r and s
    endfor
endfor

# NOTE: from the last loop, s = 36 and r = 19
# write value added, employment and GDP for years 1 and 0
varname = "v0";
var = reshape(v0,36,65);
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

varname = "v1";
var = reshape(v1,36,65);
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

varname = "e0";
var = reshape(e0,36,65);
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

varname = "e1";
var = reshape(e1,36,65);
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [var(:,s),var(:,r)]', ",", "-append")

varname = "y0";
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [y00000(s);y00000(r)], ",", "-append")

varname = "y1";
dlmwrite (filename, varname, ",", "-append")
dlmwrite (filename, [y11111(s);y11111(r)], ",", "-append")

# clear all data
clear
# END of script
